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ABSTRACT 


During the shutdown of the space shuttle main engine, oxygen 
flow is shut off from the fuel prebumer and helium is used to push 
the residual oxygen (between the oxidizer valve and the prebumer 
combustion chamber) into the combustion chamber. During this process 
a low frequency combustion instability, or chug, occurs. This chug 
has resulted in damage to the engine's augmented spark igniter due 
to backflow of the contents of the preburner combustion chamber into 
the oxidizer feed system. 

To determine possible causes and fixes for the chug, the fuel 
prebumer was modelled as a heterogeneous stirred tank combustion 
chamber, a variable mass flow rate oxidizer feed system, a constant 
mass flow rate fuel feed system and an exit turbine. Within the 
combustion chamber gases were assumed perfectly mixed. To account 
for liquid in the combustion chamber, a uniform droplet distribution 
was assumed to exist in the chamber, with the mean droplet diameter 
determined from an empirical relation. A computer program was 
written to integrate the resulting differential equations. 

Because chamber contents were assumed perfectly mixed, the fuel 
preburner model erroneously predicted that combustion would not take 
place during shutdown. The combustion rate model was modified to 
assume that all liquid oxygen that vaporized instantaneously 
combusted with fuel. Using this combustion model, the effect of 
engine parameters on chamber pressure oscillations during the SSME 
shutdown was calculated. These studies showed that decreasing the 


pressure downstream of the prebumer's exit turbine, decreasing the 
fuel temperature, increasing the helium 

temperature and decreasing the length of the line connecting the 
helium storage tank to the helium check valve all had adverse 
effects on engine stability. 
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CHAPTER I 


INTRODUCTION 

For start-up and normal operation, the Space Shuttle Main 
Engine (SSME) fuel prebumer maintains stable operation. During 
shutdown, however, both in ground test firings and actual use, the 
fuel and oxidizer prebumer s frequently undergo low frequency 
combustion instability, or chug. The chug is manifest as pressure 
oscillations in the prebumers and their feed systems. These 
oscillations have resulted in undesirable turbopump bearing loads 
find have damaged the augmented spark igniter oxidizer supply piping 
due to backflow of fuel and oxidizer from the prebumer combustion 
chamber into the feed system. This study is concerned with modelling 
the SSME fuel preburner and predicting its behavior during chug. 

Following is an overview of the SSME and fuel prebumer, an 
account of the engine shutdown process and a description of the 
model to be used. 

1.1 THE SSME AND FUEL PREBURNER 

Figure 1.1 shows a schematic of the SSME. Liquid hydrogen, 
stored at 19 K (34°R) and 0.21 MPa (30 psia), and liquid oxygen, 
stored at 90 K (164°R) and 0.69 MPa (100 psia), are the engine's 
propellants. Both the hydrogen and the oxygen enter the engine via 
low pressure pumps. These raise the pressure of the propellants to 
prevent cavitation in the subsequent pumps. After leaving the low 
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Figure 1.1: Space Shuttle Main Engine Propellant Flow Schematic 

(From Sutton, 1986) 



pressure pumps, the hydrogen and oxygen are pumped to high pressures 
by the turbopumps driven by the exit flow of the fuel and oxidizer 
preburners. At full power level the hydrogen leaves the fuel 
turbopump at 48.25 MPa (7000 psia) and the oxidizer leaves the 
oxidizer turbopump at 30.33 MPa (4400 psia). 

Most of the oxygen (86% at full power level) flows directly to 
the main combustion chamber. The remainder flows to the two 
preburners. Upon leaving the high pressure pump, 76% (at full power 
level) of the hydrogen flows directly to the prebumers. The 
remainder flows through the engine, cooling the preburner combustion 
chambers, the main combustion chamber and the nozzle. 

The oxidizer and fuel preburners serve two purposes in the 
SSME: to drive the turbines which supply power to the high pressure 
pumps and to preheat the hydrogen for combustion in the main 
combustion chamber. At full power level, the fuel preburner operates 
at around 1090 K (1960°R) and 39.01 MPa (5660 psia) and the oxidizer 
preburner operates at around 830 K (1500°R) and 38.60 MPa (5600 
psia). During the fuel preburner chug, the pressure in the fuel 
preburner varies around 4.48 MPa (650 psia). Actual average chamber 
temperature data is not available, but the Marshall Space Flight 
Center SSME transient model predicts mean temperatures between 440 K 
and 800 K (800°R and 1400°R) during the chug. 

Injection of liquid oxygen and supercritical hydrogen into the 
fuel prebumer is accomplished by means of 264 coaxial injectors. To 
cool the plate, hydrogen also is injected through small orifices 


* 


m 


across the face plate. At full power level, just before injection, 
the hydrogen is at about 160 K (284°R) and 43.08 MPa (6250 psia) and 
the oxygen is at about 120 K (214°R) and 48.4 MPa (7020 psia). 

At full power level, the fuel preburner operates at an 
equivalence ratio ((fuel/oxidizer)/(fuel/oxidizer) sto ^ c k£ 0metr £ C ) of 
about 8. Mass flow rates of hydrogen and oxygen into the fuel 
preburner are 36.4 kg/s (80 lb/s) and 37.7 kg/s (83 lb/s), & 

respectively. Hot gases from the preburners exit the preburner 
combustion chambers through turbines and enters the main combustion 
chamber via the hot gas manifold and the main injection plate. In 
the main combustion chamber the products from the two preburners are 
burned with oxygen. The equivalence ratio in the main combustion 
chamber is about 1.33. 


1.2 SSME SHUTDOWN AND THE FUEL PREBURNER CHUG 

When the SSME shuts down, it is important that combustion 
ceases in the preburners and that oxidation does not begin in the 
feed system, so a helium purge is used to push oxygen from the feed 
system and preburners and assure that the system is above the rich 
limits of combustion of hydrogen with oxygen. It is during the 
helium purge that the fuel prebumer chug takes place. Figure 1.2 
shows the fuel prebumer manifold. The helium purge orifice is a 
small pipe that links the helium supply piping to the pipes that 
reach the helium check valves. The helium check valves restrain the 
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helium until the pressure downstream of the valves drops below 5.17 
MPa (750 psia). 

Shutdown is accomplished and controlled by the closing of the 
oxidizer prebumer oxidizer valve, the fuel prebumer oxidizer 
valve, the main oxidizer valve and the main fuel valve and the 
opening of the helium purge check valves. Table 1.1 (adapted from 
George, 1984) gives the valve sequence and the resulting events 
during the SSME shutdown. Two points should be noted from Table 1.1. 
First, before chugging begins in the fuelside preburner, both the 
fuelside preburner oxidizer valve and the main oxidizer valve are 
closed; the fuelside preburner is isolated from the oxidizer feed 
system. Second, no valve openings or closings coincide with the 
onset or ending of chug. These points are further discussed in 
Chapter II. 

When the fuelside preburner chug begins is dependent on the 
power level from which the engine is cut off and the diameter of the 
helium purge orifice. Pressure traces for typical fuelside preburner 
chugs for SSME test firings are found in Figure 1.3. These traces 
are plots of the deviation of the chamber pressure from the mean 
chamber pressure versus time. The traces were filtered to show only 
frequencies between 50 and 200 Hz. Figure 1.3 (a) is a pressure 
trace for a small helium purge orifice diameter (0.17 cm [0.068 
inches]) and Figure 1.3 (b) is a trace for a large helium purge 
orifice diameter (0.74 cm [0.291 inches]). Using the small helium 
purge orifice resulted in a short chug whose amplitude is about 2.89 


Table 1.1: Sequence of Events During the SSME Shutdown 


TIME AFTER 
CUTOFF (SEC) 


1.4 


1.8 


2.2 


2.3 


2.3 - 2.5 


2.5 


3.0 


EVENT 


OXIDIZER PREBURNER OXIDIZER VALVE SEALS CLOSE 

OXIDIZER PREBURNER PRESSURE DROPS SMOOTHLY FROM 
800 PSIA TO 500 PSIA 

HELIUM CHECK VALVES OPEN, PUSHING OXIDIZER 
FROM THE OXIDIZER PREBURNER MANIFOLD INTO 
COMBUSTION CHAMBER 

OXIDIZER PREBURNER CHAMBER PRESSURE RISES TO 
750 PSIA AND STABILIZES. 

HELIUM FLOW SUBSIDES AFTER THE INITIAL SPIKE 

FUEL PREBURNER OXIDIZER VALVE SEALS CLOSE 

HELIUM FLOW RATE IMMEDIATELY RISES 

MAIN OXIDIZER VALVE CLOSES (AFTER THE FUEL 
PREBURNER OXIDIZER VALVE) 

HELIUM FLOW STABILIZES, AS DO THE FUEL PREBUR- 
NER AND OXIDIZER PREBURNER CHAMBER PRESSURES 
FOR THE SMALL ORIFICE TESTS 

HELIUM FLOW RATE CONTINUES TO INCREASE FOR THE 
LARGE ORIFICE TESTS 

CHUGGING BEGINS IN THE FUELSIDE PREBURNER; 

LARGE AMPLITUDE FOR THE LARGE ORIFICE 
AND SMALL AMPLITUDE FOR THE SMALL ORIFICE 

LOW LEVEL CHUGGING MAY ALSO BEGIN IN THE OXI- 
DIZER PREBURNER 

LARGE AMPLITUDE CHUG BEGINS FOR THE SMALL 
HELIUM PURGE ORIFICE TESTS 

A TEMPORARY REDUCTION IN THE AMPLITUDE OF THE 
THE FUELSIDE PREBURNER CHUG OCCURS DURING 
THE LARGE ORIFICE TESTS 

CHUGGING ENDS 


3.7 


THE MAIN FUEL VALVE CLOSES 





DEVIATION FROM THE MEAN 
DEVIATION FROM THE MEAN CHAMBER PRESSURE (psi) 

CHAMBER PRESSURE (psi) , 


» 


I 




(b) Helium Purge Orifice = 0.291" 

Figure 1.3: Pressure Traces for the SSME Fuel Preburner Chug 

(From George, 1984) 

(Pressure Traces Were Filtered to Show Only Frequencies 
Between 50 and 200 Hz) 
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MPa (420 psi). Using the large orifice resulted in a longer but less 
severe chug. About three seconds after cutoff the chug ends and the 
chamber pressure and temperature decrease smoothly. The main fuel 
valve closes after the purge is complete. 

Because the chug takes place after the fuel preburner and its 
piping are isolated from the oxidizer feed system and before the 
main fuel valve is closed, George (1984) concluded that the chug was 
not directly related to or triggered by valve openings and closings 
and that it was associated with the oxidizer purge process. As 
helium replaces oxygen in the oxidizer feed system, the fluid in the 
pipes becomes more compressible and propellant feed rates to the 
preburners become more sensitive to chamber pressure. This effect is 
described in greater detail in Chapter II. 

1.3 MODELLING THE SSME FUEL PREBURNER AND FEED SYSTEM 


Because of high mass flow rates, very reactive propellants, 
high chamber pressure and the fact that combustion takes place on 
many diffusion flames surrounding droplets throughout the combustion 
chamber, combustion in a liquid propellant rocket engine can be 
considered to take place homogeneously through the combustion 
chamber (Summerf ield, 1951). Consequently, the present study 
modelled the fuelside prebumer combustion chamber as a stirred tank 
reactor; contents are considered well mixed and properties are 
considered uniform throughout the chamber. Liquid droplets were 
distributed evenly through the combustion chamber and all had the 
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same characteristic diameter. The heterogeneous stirred tank reactor 
is more fully described in Chapter III. 

The governing equations for the combustion chamber were derived 
from the conservation of mass applied to the liquid in the 
combustion chamber, the conservation of mass applied to the gases in 
the combustion chamber, the conservation of species, the 
conservation of energy and the perfect gas law. Auxiliary relations 
used to obtain terms in the governing equations included empirical 
expressions for average droplet diameter, droplet evaporation rate, 
thermodynamic properties for the gases in the combustion chamber and 
elementary reaction rates. Exhaust mass flow rate from the preburner 
was given as a function of the ratio of the pressure downstream of 
the exit turbine to the combustion chamber pressure. At each time, 
combustion chamber pressure, temperature, density (excluding liquid 
mass, which is assumed to occupy negligible volume), liquid density 
(kg of liquid per unit volume of the combustion chamber) and mole 
number of each species (kgmoles of the species per total mass (kg) 
in the combustion chamber) were calculated. 

In the feed system the effects of the compressibility of helium 
were important, as was the feed system geometry, so a multiple pipe, 
multiple node feed system containing distinct liquid and gaseous 
phases was used to predict the transient feed system behavior. 
Governing equations for the feed system included the conservation of 
momentum along each pipe and the conservation of mass at each node. 
Time dependent variables were the fluid velocity in each pipe, the 
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liquid/ vapor interface position in each pipe and the density at each 
node. 

The combustion chamber and feed system governing equations were 
solved numerically on a digital computer, integrating over time. 
Since the equations are stiff (i.e., characterized by widely 
differing time constants for~the dependent variables (Pratt and 
Radhakrishnan, 1986)) steps were taken to make the equations 
solvable without using an excess of computer time. In addition to 
using an integration routine written for solving systems of stiff 
differential equations, a steady state approximation was made for 
the concentrations of radicals and atoms in the combustion chamber 
to facilitate the solution of the governing equations. 

Derivation and solution of the governing equations is given in 
Chapter III. 

Success of this model is determined by its ability to predict a 
chug as encountered in the SSME fuel preburner during shutdown and 
its ability to vary engine parameters and find a way in which the 
fuel preburner chug can be lessened in severity or eliminated. 
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CHAPTER II 


LITERATURE SURVEY 

2.1 OVERVIEW OF DIFFERENT INSTANCES OF COMBUSTION INSTABILITY, THEIR 
CAUSES AND EFFECTS. 

The phenomenon of combustion instability exists to some degree 
in all combustion processes . The sputtering of a match, the 
crackling of a log on a fire and the flickering of a candle are all 
everyday examples of combustion instability. Such early scientists 
as Rayleigh and Farraday noted pressure oscillations as sound waves 
generated during the burning of substances in a tube. Rayleigh 
(1945) described these oscillations as 

... vibrations maintained by heat, the heat being communicated 
to the mass of air confined in the sounding tube at a place 
where, in the course of vibration, the pressure varies. 

He concluded that "In consequence of the variable pressure within 
the resonator, the issue of gas, and therefore, the development of 
heat, varies during vibration." This is an apt description of the 
processes involved in a periodic combustion instability. 

Combustion instability can be classified in a number of 
different ways. First, it can be either periodic or random. An 
example of a periodic instability is a pressure wave travelling 
around a combustion chamber due to varying combustion rates. Spikes 
or pops (spontaneous explosions) in rocket combustion chambers are 
examples of random instability. Second, combustion instabilities are 


either chamber, system or intrinsic instabilities (Williams, 1985). 
Chamber instabilities involve only the processes that occur within a 
combustion chamber, while system instabilities involve interaction 
between processes occurring in a combustion chamber with processes 
occurring elsewhere in the system, such as in propellant feed lines 
or at the combustion chamber exit. An intrinsic instability involves 
only the combustion process and occurs whether or not there is a 
combustion chamber. Finally, periodic combustion instability can be 
classified according to its frequency. In rocket motor combustion 
instability, there exist three fairly distinct instability frequency 
ranges. High frequency instabilities (normally greater than 1000 Hz) 
are chamber instabilities involving longitudinal, transverse or 
radial waves within the rocket combustion chamber. Low frequency 
instabilities (normally less than 300 Hz) are system instabilities 
which couple changes in the propellant feed rate with changes in 
combustion rate in the combustion chamber. Intermediate frequency 
instabilities (frequencies between 300 Hz and 1000 Hz) involve 
aspects of both system and chamber instabilities: waves propagate 
within the combustion chamber and also within the feed system. 

Periodic combustion instability is observed in ramjet engines 
and solid rocket engines as well as in monopropellant and 
bipropellant liquid rocket engines. Low frequency combustion 
instability in ramjet engines is a result of vortices in the 
combustion chamber and is very dependent upon chamber geometry (Yang 
and Culick, 1986). There are a number of parallels between liquid 
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rocket and solid rocket engine low frequency combustion 
instabilities. Schoyer (1986) describes the physical processes 
involved in solid rocket chugging as follows: 


The basic idea is that due to oscillatory combustion there is a 
fluctuating heat transfer into the propellant which 
superimposes an oscillation on the steady state profile in the 
pyrolizing propellant. Because the solid propellant depends on 
the (fluctuating) surface temperature, the pyrolysis rate will 
also vary, causing a fluctuating mass flow into the combustion 
chamber. If this fluctuating mass flow is more or less in phase 
with the pressure fluctuations, the two effects may enhance 
each other (causing resonance) and lead to oscillatory 
combustion. 


Just as heat transfer and pyrolysis provide gaseous reactants 
in the burning of solid propellants, heat transfer and vaporization 
provide the gaseous reactants in the burning of liquid propellants. 
Combustion in solid and liquid rockets differs in two important 
ways. First, combustion can be considered to take place on a single 
planar flame front in a solid rocket combustion chamber, whereas 
combustion in a liquid propellant combustion chamber involves many 
diffusion flames that surround droplets. Second, the liquid 
propellant feed rates must also be considered for a liquid rocket 
engine. 

As mentioned above, liquid rocket engines can undergo low, 
intermediate and high frequency combustion instability. The physical 
manifestations of these include noise, pressure and temperature 
oscillations in the combustion chamber (longitudinal, radial and 
transverse), changes or even reversals in the mass flow rate of the 
feed system, increased rates of heat transfer to engine components. 



decreased engine performance and possibly severe vibrational loads 
on engine components. 

For a liquid rocket engine, the most destructive type of 
combustion instability is high frequency instability, sometimes 
called screech, scream or acoustic instability. Screech can result 
in vibrational loads on rocket components and, more importantly, 
increased rates of heat transfer to chamber walls. Lawhead and Combs 
(1963) observed burnout of the injector face of a liquid rocket 
combustion chamber between 200 and 300 msec after the beginning of 
screech. Penner and Datner (1955) stated that performance might be 
improved by screech, although this would not offset the negative 
effects of screech. 

During screech oscillations of pressure and temperature usually 
take place only within the combustion chamber of a rocket motor, 
although oscillations can also occur in the feed system. 

(Coultas, 1972) . Upon a perturbation, thermodynamic properties and 
chemical composition can change in a part of the combustion chamber, 
leading to a local change in the rate of combustion (Crocco et al., 
1960). This change initiates a wave which propagates through the 
combustion chamber. The wave can be linear (discrete) or nonlinear 
(coupled with other waves) and radial, longitudinal or 
circumf erential. The propellant feed system is usually not directly 
involved in the propagation of pressure and temperature waves, but 
is important in that it determines the amount and distribution of 
reactants in the combustion chamber. The mass flow rate in the feed 
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system does not oscillate during screech because oscillations in the 
chamber take place at a higher frequency than those to which the 
feed system can respond. 

Among possible sustaining mechanisms for screech are the 
physical delays associated with the combustion process (atomization, 
vaporization and reaction rate), the sensitivity of the rate of 
combustion to pressure and temperature changes and the "explosion" 
of liquid droplets which are heated to temperatures above their 
critical temperatures (Coultas, 1972). Reardon, McBride and Smith 
(1966) showed that spatial nonuniformities in propellant injection 
rates across the injector were also, at times, responsible for 
screech. Their work was verified by experiment. 

To prevent destructive high frequency combustion instability, 
liquid propellant rockets may need baffles or liners to cause 
combustion chamber conditions (pressure, temperature, distribution 
of propellants) to be more uniform or to change the resonant 
frequency of the engine (Williams, 1985). A change in the propellant 
spray field can also prevent destructive scream (Coultas, 1972). 

Both of these fixes involve local control of combustion rates. 

Less common in liquid rocket motors, intermediate frequency 
instability, or buzz, involves waves in the propellant feed system 
as well as spatial combustion chamber pressure variations (Coultas, 
1972). Buzz differs from screech in that it involves lower frequency 
disturbances and oscillations in the feed system. It differs from 
low frequency instability in that spatial variations in pressure. 
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temperature and chemical composition exist in the combustion 
chamber. Buzz causes structural loads on engine components and 
standing waves in the feed system pipes and may initiate high 
frequency instability. 

Scream involves pressure and temperature waves that propagate 
in the combustion chamber as a result of spatially varying rates of 
combustion. Buzz is characterized by standing waves in the 
propellant feed system and spatial variations in combustion rate in 
the combustion chamber. In contrast to scream and buzz, combustion 
chamber pressure can be considered uniform during low frequency 
instability, or chug. In chug, aggregate oscillations in chamber 
pressure couple with oscillations in propellant feed rate, resulting 
in loss of performance, structural loads and possible failure of 
engine parts. The coupling of oscillations in chamber pressure and 
propellant feed rate is a result of the time delay between the 
injection of propellants into the combustion chamber and their 
combustion. Chugging in liquid bipropellant rocket engines is 
further described in section 2.2. 

2.2 LOW FREQUENCY COMBUSTION INSTABILITY IN LIQUID BIPROPELLANT 
ROCKET ENGINES 

Chugging in liquid bipropellant rocket engines is characterized 
by pulsations in the rate of propellant injection into the 
combustion chamber and variations in combustion chamber pressure 
(Penner and Datner, 1955). These pulsations are generally taken to 
be the result of the time lag that exists between the injection of 


liquid propellants into the combustor and their conversion to hot 
product gases, as proposed by Summerf ield (1951). Because of the 
inertia of the propellant mass in the feed lines, mass flow rates 
into a combustion chamber are not immediately responsive to changes 
in chamber pressure. Because of the combustion time lag, combustion 
chamber temperature, pressure and composition are not immediately 
responsive to changes in propellant mass flow. When response delays 
in the feed system and combustion system become coupled, large 
amplitude chugs can result. Figure 2.1 shows the coupling of the 
feed system and combustion chamber for a chugging rocket engine. 

Because the conversion of low temperature propellants to high 
temperature products in a liquid bipropellant rocket engine is such 
a complex process, it would be difficult for a chug model to 
incorporate all facets of the process. An effort must be made to 
determine the most important parts of the process. In the SSME 
fuelside preburner, liquid oxygen and gaseous hydrogen combust. This 
involves injection of the hydrogen and oxygen, atomization and 
vaporization of the oxygen, mixing of the gaseous hydrogen and 
oxygen vapor and chemical reaction of the gaseous species. In his 
study of screech in liquid bipropellant rockets, Priem (1966) took 
atomization, vaporization and chemical reaction as the most 
important processes in combustion instability. In the current study, 
the author will also take these three processes to be the most 
important parts of the combustion process and neglect other parts 
such as gas phase mixing and wall catalyzed reactions. Atomization, 
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vaporization and chemical reaction and their effect on chugging are 
briefly described below. 

The atomization of a propellant stream is primarily dependent 
upon the way the stream is injected into the combustion chamber. 
Propellants can be premixed or distinct upon injection and injection 
can be accomplished by parallel jets, impinging streams or coaxial 
streams. After injection, streams may be deflected by a splashplate. 
The SSME fuelside prebumer employs coaxial injectors. Coaxial 
injectors are primarily used when one of the injected fluids is 
gaseous and the other is a liquid, since they are well adapted to 
mixing a gaseous stream with a liquid one (Dykema, 1972). 

During atomization, streams break into ligaments which break 
into droplets. Droplets can either evaporate, shatter or interact 
with each other. Priem (1966) said that atomization rate is 
proportional to the mass of unatomized liquid in the combustion 
chamber and a less strong function of chamber density and axial 
velocity of the liquid propellant upon injection. The mass of 
unatomized propellant in the combustion chamber at a given instant 
is a function of the injection rate previous to that time, so a 
change in injection rate of propellant will affect the atomization 
rate only after a time delay. 

Spatial injection droplet distribution has been shown to have a 
great effect on screech (Reardon et al., 1966), but does not greatly 
affect chugging, since it is aggregate chamber pressure oscillations 
and not spatial pressure oscillations that are able to couple with 


the feed system. Webber (1972) showed that gross changes (for 
example, not including any small droplets in the droplet 
distribution in the combustion chamber) in the estimation of droplet 
sizes in a chug model drastically changed his predicted pressure 
oscillation amplitudes and frequencies. 

The vaporization rate for a liquid propellant in a rocket 
engine combustion chamber is dependent on the atomization process 
(droplet size and surface area), the heat transfer rate to the 
droplet and the rate at which propellant diffuses away from the 
droplet. To prevent chug, the part of the combustion time delay 
associated with vaporization should be made as short as possible. In 
the SSME fuels ide preburner, coaxial injectors inject high speed 
hydrogen around a jet of liquid oxygen. The greater the speed of the 
hydrogen, the more shear on the liquid oxygen jet, the finer the 
atomization, the greater the vaporization rate and the shorter the 
vaporization delay time. Higher chamber temperatures also result in 
greater rates of vaporization of liquid propellants due to greater 
rates of heat transfer to droplets. Another way to increase the 
vaporization rate is to increase the temperature of the liquid 
propellant. For higher temperature liquid propellants, less energy 
must be transferred to the droplets to raise their temperature and 
vaporize them. 

Chemical reaction rates, or equivalently, propellant 
combinations, can greatly influence combustion stability in a rocket 
engine, though often chemical reaction rate is taken as a less 


important factor (Priem, 1966, and Reardon, 1972). Reaction rate can 
be important if it is sufficiently slow and makes up a significant 
fraction of the combustion time delay. Summerfield (1951) suggested 
that chug could be eliminated in some engines through the use of 
either more reactive propellants or some sort of catalysis. 

A number of experimental investigations have been conducted to 
determine how rocket design parameters such as combustion chamber 
pressure, propellant combination, injector pressure drop, chamber 
length, number of injectors, feed system mass flow rates and 
propellant mixture ratio influence chugging. The results of two 
papers will be highlighted below. 

Both Barrere and Moutet (1956) and Heidmann et al. (1967) noted 
that higher chamber pressure leads to higher frequency but lower 
amplitude pressure oscillations (Figures 2.2(a) and 2.2(b)). 
Temperature and pressure oscillations were roughly in phase during 
chug (Barrere and Moutet). 

As characteristic chamber length (usually defined as the 
chamber volume divided by the area of the exit nozzle throat) is 
increased, the range of stable operating conditions is decreased 
(Heidmann et al.) while the frequency of pressure oscillations 
decreases (Barrere and Moutet). Chambers of the same characteristic 
length but different volumes were studied by Barrere and Moutet. 

Both chambers exhibited roughly the same pressure oscillation 
frequency, with the smaller chamber undergoing oscillations of 
slightly higher amplitudes. 



(a) Pressure Oscillations for an Average 
Chamber Pressure of 0.88 MPa 



(b) Pressure Oscillations for an Average 
Chamber Pressure of 1.18 MPa 


Figure 2.2: Effect of Chamber Pressure on Chug Amplitude and 
Frequency (Barrere and Moutet, 1956) 
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Increasing the number of jets (injectors) of a combustion 
chamber decreased the range of stable operating conditions (Heidmann 
et al.)» pres uma bly because of droplet and stream interactions 
driven by the pressure changes. 

Changing the mixture ratio within a moderate range (0.9 < m/m s 
< 1.4 , where m is the fuel/oxidizer mixture ratio and m s is the 
stoichiometric mixture ratio) had very little effect on chug 
frequencies for an engine whose propellants were nitric acid and 
furfurylic alcohol. The lowest amplitude pressure oscillations 
occurred for a stoichiometric mixture. This is due to a reduced 
combustion delay time, since collisions between appropriate 
molecules of reactants are most likely for a stoichiometric mixture. 

2.3 CONTROLLING LOW FREQUENCY COMBUSTION INSTABILITY 

Having described the processes involved in chugging, a number 
of methods for controlling chugging will be presented. First, 
combustion chamber geometry (characteristic length) can be changed. 
In many cases, as with the SSME, chug is encountered only after the 
engine has been built, so changing the chamber geometry is not a 
viable alternative. In addition, since chugging only occurs in the 
SSME during shutdown, extensive changes in chamber geometry is not 
merited. 

Second, the pressure drop across the injectors can be 
increased. A larger pressure drop across the injectors means that a 
change in downstream (chamber) pressure will be a smaller fraction 


of the upstream pressure (behind the injector) and will not effect 
feed system mass flow rates very much. Friction in the feed system 
and injectors damps oscillations in propellant flow. Any increase in 
the energy of the propellant stream will be reduced by the feed 
system and injector friction. Consequently, choosing an injector 
with much friction could result in a more stable engine. On the 
other hand, a greater pressure loss for the stream entering the 
combustion chamber would be the cost. 

2.4 THE SSME FUELSIDE PREBURNER CHUG 

There are a number of factors peculiar to the SSME fuelside 
preburner that contribute to the onset of chugging. First, as 
described in the introduction, the chug occurs during the SSME shut- 
down while helium progressively replaces liquid oxygen in the oxi- 
dizer feed system. The presence of helium in the oxidizer feed sys- 
tem makes the oxidizer mass flow rate more sensitive to changes in 
the chamber pressure than if there were no helium in the feed sys- 
tem. Gaseous helium is more compressible than liquid oxygen, so an 
increase in chamber pressure results in compression of the helium 
and a greater reduction in the oxidizer mass flow rate than if no 
helium were present. For a temperature and pressure near those at 

which the chug occurs (5.0 MPa and 100 K), the density and bulk 

3 3 

modulus of oxygen are about 1000 kg/m (62.4 lbm/ft ) and 221 MPa 

2 

(32,000 psia), respectively. For helium they are about 17 kg/m 

3 

(1.06 lbm/ft ) and 12.4 MPa (1800 psia), respectively. 


In the SSME fuelside preburner, liquid oxygen is burned at a 
low equivalence ratio (about 8.0 at rated power level) with 
hydrogen. Hydrogen has a very low critical temperature (33.3 K 
[60°R]) and pressure (1.29 MPa [188 psia]), so it enters the 
combustion chamber as a supercritical gas. The critical pressure of 
oxygen is about 5.09 MPa (738 psia), so during the SSME fuelside 
preburner chug, which takes place at pressures which vary around 
4.83 MPa (700 psia), oxygen is a near critical liquid. 

The conversion of gaseous hydrogen and oxygen to water vapor 
occurs as a chain reaction involving the radicals and atoms OH, 0, 

H, HO 2 and ^ 2^2 (Kuo, 1986). At high pressures, the radicals HO 2 and 
H 2 O 2 dissociate and for high temperatures and high mass flow rates, 
reactions at the wall are not very important. Since the SSME 
fuelside preburner operates at high pressures, mass flow rates and 
temperatures, the following reactions can be considered to be the 
elementary steps in the overall conversion of and O 2 to water 
vapor : 


• 

H 2 + 0 2 

--> 

2 OH 



H 2 + M 

--> 

2H ■ 

1- M 


OH + H 2 

--> 

h 2 o 

+ H 

• 

H + 0 2 

--> 

OH 

+ 0 


0 + h 2 

--> 

OH 

+ H 


H + H + M 

--> 

»2 

+ M 

• 

H + 0 + M 

--> 

OH 

+ M 


0 + 0 + M 

--> 

°2 

+ M 
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H + OH + M --> H 2 0 + M 

where M represents a third body, inert to the reaction . For a 
rocket engine combustion chamber, which operates at high 
temperatures and pressures, the atoms and radicals produced in the 
chain branching reactions will be consumed in the termination 
reaction nearly as soon as they are produced. 

George (1984) observed during test firings of the SSME that 
using a 0.068 inch diameter helium purge orifice for the SSME 
fuelside preburner as opposed to a 0.291 inch orifice resulted in a 
34% increase in maximum peak to peak amplitudes of pressure 
oscillations (from 335 psia) and a 38% decrease in chug duration 
(from 0.72 sec). Changing the purge orifice diameter did not 
substantially change the chug's center counted frequency (the 
average frequency at the chug's temporal center). The center counted 
frequency for the small and large orifices were 117 Hz and 117.5 Hz, 
respectively. Traces of pressure for small and large orifice engine 
test firings are shown in Figures 1.3(a) and 1.3(b) (page 8). Chug 
ending time was not affected by the change of the helium purge 
orifice diameter, but chug starting time was. Because of this, 

George ruled out oxygen depletion as the chug terminating mechanism. 
Helium flow rate and compressibility are evidently important in 
determining the character of the SSME fuelside prebumer chug. 

The only opening or closing of a valve that occurs near the 
chug's end is the closing of the main oxidizer valve. As described 
in Chapter I, during chug the fuelside preburner is isolated from 
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the oxidizer feed system, so valve openings and closings are 
probably not a concern in the SSME fuels ide preburner chug. The 
valves that restrain the helium before the purge begins are poppet 
type valves that remain open even when the pressure drop across the 
valve is less than the cracking pressure (George, 1985). Therefore, 
helium purge valve clatter probably does contribute to the chug. 

Because of the fuelside preburner turbine blade inertia, the 
exit mass flow rate does not respond immediately to a change in 
combustion chamber pressure. A relation for exit mass flow rate as a 
function of chamber pressure and temperature and pressure downstream 
of the turbine was given by Nguyen (1981) and is found in Chapter 
III. 

Because of the above considerations, the model of the SSME 
fuelside preburner allows for the compressibility of helium in the 
piping system, a combustion time lag due to the atomization and 
vaporization times for oxygen and turbine blade inertia. 

2.5 ROCKET ENGINE STABILITY PREDICTION MODELS 

Most chug models have attempted to predict the domain of 
operating conditions within which a rocket engine cam operate 
without undergoing destructive oscillations in pressure and 
temperature. By varying inputs such as combustion chamber pressure, 
injector pressure drop, combustion time delay or instability 
frequency to these models, stability boundaries for specific engines 
and general criteria for rocket engine stability can be found. 


Figure 2.3 shows the rocket engine model used by Summerfield 
(1951). He related the rate of change of combustion chamber pressure 
with respect to time, t, at a given instant, to the propellant 
injection rate at a time t- 8 , using the equation: 


dP 
c 

dt 


R T 

* -^PcV2 (t ' e) 


R T 

c 

L*c* 


P 

c 


[ 2 . 2 ] 


where T c is the combustion chamber temperature (K),p c is the 

3 

combustion chamber density (kg/m ), V is the combustion chamber 
3 

volume (m ), L* is the characteristic chamber length (m), c* is the 
characteristic propellant velocity (m/s), R is the gas constant for 
the gases in the combustion chamber (J/kg*K), A 2 is the area (m 2 ) at 
station 2 in Figure 2.3 and V 2 (t- 8 ) is the propellant velocity (m/s) 
at station 2 at time t-0. This implies that a mass of propellant 
injected into the combustion chamber occupies negligible volume and 
does not react at all until 8 seconds later when it is converted to 
products instantaneously. 8 is the combustion time delay and can be 
a function of chamber pressure, temperature and composition and of 
propellant composition and velocity. Characteristic chamber length, 
L*, is defined as the volume of the chamber divided by the average 
cross sectional area. Summerfield defined the characteristic 
propellant velocity, c*, as 

c* = f(y)-^ R t c 

where Y is the ratio of the specific heats in the chamber . 
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Figure: 2.3 Schematic of Summerf ield's Liquid Propellant 
Rocket System (Summerf ield, 1951) 
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Summerfield used conservation of energy for the control volume 
bounded by surfaces and along with continuity and equation 
[2.2]. Differentiating the result with respect to time and 
linearizing by saying that V 2 (t) - v^Ct-©) is small, gave an 
equation of the form 

u ,, ( t) + A u'(t) + B u(t) + C u(t-e) = 0 [2.3] 


where u = d(v 2 )/dt. Solutions of [2.3] are of the form 

u = £ U n exp((A n + 6> n )t) [2.4] 

Substituting [2.4] into [2.3], the values of A n and which satisfy 
[2.3] can be found. 

Analyzing the eigenvalues of [2.3], Summerfield determined that 
the following condition must be met for a liquid rocket engine to be 
stable (i.e., to not undergo undamped pressure oscillations): 



P 

c 



+ 


c*L* 2 AP 



e 


[2.5] 


where AP is the pressure difference between the propellant storage 
tank and the combustion chamber (psi) and is the length of the 
propellant feed line (Figure 2.3). 

Inequality [2.5] indicates that an unstable rocket engine can 
be made stable by 

1. increasing the feed system pipes length, 

2. increasing the propellant feed mass flow rate. 
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3. increasing the pressure drop from the storage tank to the 
combustion chamber, 

4. decreasing the combustion chamber pressure or 

5. increasing the characteristic chamber length. 

Evaluation of [2.5} with the appropriate SSME fuels ide 

preburner shutdown values (1 2 = 0.5 m [1.64 ft], m = 33.91 kg/s 
[74.8 lb/s], P c « 4.825 MPa [700 psia], A 2 = 0.0025 m 2 [0.027 ft 2 ], 
c* = 60 m/s, [197 ft/s] [estimated from calculations by VanOverbeke 
and Claus, 1986], L* - 0.1 m [0.328 ft], R = 2080 J/kg-K [0.497 
Btu/lbm* °R] , AP = 0.345 MPa [50 psia] and T c = 700 K [1260°R]) 
yields the stability criterion 

9 < 0.0014 s 

(i.e., the combustion time delay must be less than 440 ps). 

Crocco (1951) began his stability analysis for a monopropellant 
rocket engine using conservation of mass in the combustion chamber. 
Assuming, like Summerfield, that propellants have negligible volume 
until burned and that the mass burned over the time interval dt is 
equal to the mass of propellant injected over the time interval 
d(t-Q), continuity for the combustion chamber can be written 

dM d9 

— ^ + m (t) = (1 - — ) m.(t-8) [2.6] 

dt ® dt 1 
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where m and m. are the mass flow rates out of and into the 
e 1 

combustion chamber and M is the mass of product gases in the 

s 

combustion chamber. 

Crocco divided the combustion delay time, 8, into an 

insensitive part, 9 , which does not vary with small variations in 

8 

chamber conditions and a sensitive part, z . 


0 = 8 + z 

g 


[2.7] 


9 can be given, for instance, by L*c*/R*T . Crocco assumed that z 

g c 

obeys the relation 


z P n = const. [2.8] 

c 

Using [2.8] and [2.7], Crocco nondimensionalized and linearized 

[ 2 . 6 ]. 

Applying the conservation of mass and momentum to the feed 
system, piping system behavior was seen to be dependent on three 
dimensionless groups: one involving the pressure drop across the 
injectors, the second involving the inertia of the propellant in the 
pipes and the third giving the relative elasticity of the pipes and 
the propellant within. 

Using a version of [2.6] and the feed system equations, Crocco 
formulated a characteristic equation for a monopropellant rocket 
engine and found the stability boundaries that corresponded to 
various values of n in equation [2.8]. He concluded that using a 
varying time lag changed stability predictions, predicting a greater 
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range of stable operating conditions, but the changes were not 
great. 

Using an annular section of a liquid bipropellant rocket engine 
combustion chamber, Priem determined the effect that atomization, 
vaporization and chemical kinetics had on screech. Though this model 
was for screech, the way that the combustion rate was modelled can 
be extended to a chug model. 

Priem calculated combustion rate in three different ways. 

First, he assumed that chemical kinetics was the most important part 
of the combustion process. Atomization and vaporization times were 
assumed short enough that the propellants that were injected into 
the combustion chamber were immediately available for chemical 
reaction. In the second model vaporization was the controlling 
process. Atomization was considered very fast, so the liquid 
injected into the combustion chamber immediately assumed a given, 
nonvarying droplet distribution profile. Chemical reaction time was 
also considered negligible. Consequently, the combustion rate was 
equal to the vaporization rate. Finally, the times required for 
atomization and chemical reaction were considered very small and the 
combustion rate was set equal to the atomization rate. 

Using conservation of mass, energy and axial momentum, Priem 
generated stability boundaries. The physical processes of 
vaporization and atomization were shown to be more important than 
chemical kinetics for screech. Stability was shown to be very 


sensitive to the mass of unburned propellant in the combustion 
chamber. Priem confirmed his stability boundaries experimentally. 

In a companion study to this one, Lim (1986) determined the 
stability boundaries of the SSME fuelside preburner during shutdown. 
Lim assumed 

1. combustion chamber pressure is uniform and oscillates about 
a steady state value, 

2. combustion chamber temperature is uniform, 

3. the combustion delay time is the same for both propellants 
and 

4. reaction rate is infinite and chemical kinetics can be 
ignored. 

Like Summerf ield, Lim assumed that the rate that propellants 
are burned is equal to their rate of injection at a previous time, 
t-8, where 


9=0 +9 [2.9] 

v m 

corresponds to vaporization time and 0^ to mixing time. 

Conservation of mass was applied to the combustion chamber, 
leading to 

— (P c V) = m fb (t-9) + m Qb (t-0) - m e (t) [2.10] 

dt 

• • 

where m^ b and m Qb are fuel and oxidizer burning rates, respectively 

(determined by conditions at time, t-9) and m is the chamber exit 

e 

mass flow at time t. The system variables m^, m Qb , m g and P c were 
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expressed in terms of a steady state value plus a perturbation and 
substituted into [2.10] which was linearized. The resulting 
characteristic equation was solved and stability boundaries were 
generated. 

In generating stability boundaries. Lim varied chamber 
pressure, fuel mass flow rate, oxidizer mass flow rate, fuel 
injection temperature and oxidizer injection temperature. Lim made 
the following conclusions from his stability boundaries. 

1. The prebumer system is generally more stable at lower 
chamber pressures. 

2. For oxidizer and fuel temperatures of 40K, the system is 
inherently unstable. For oxidizer and fuel temperatures of 
120 K, stability can be achieved only for P c < 3.79 MPa (550 
psia). Chugging could not be avoided for low fuel 
temperatures . 

3. Fuel flow rates did not influence stability. 

4. High oxidizer flow rates (or high helium purge rates) made 
the system stable. 

2.6 STIRRED TANK REACTOR MODELS 

A chemical reactor can be modelled as a batch reactor, a plug 
flow reactor or a stirred tank reactor (STR). In a batch reactor 
(Figure 2.4 (a)) reactants are fed into the reactor and left to 
react. After some time, the reactor is emptied. A plug flow reactor 
is shown in Figure 2.4 (b). Elements of fluid travel in an orderly 
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way through the chamber with no element overtaking any other element 
(i.e., no mixing in the direction of the flow, but infinite radial 
mixing). A STR is a reactor for which contents are considered well 
mixed and properties of the contents are uniform throughout the 
combustion chamber (Figure 2.4 (c)). One ramification of the STR 
assumption is that some unburned reactants will always be present in 
the chamber outflow. 

Webber (1972) used a heterogeneous STR to determine the 
transient behavior of a liquid rocket engine undergoing chug. His 
model included a simple feed system, injectors, a combustion chamber 
and a nozzle. 

A lumped parameter approach was used for the feed system, 
resulting in the equation 

^ ? t - P c - (1/2) f Pc q M [2.11] 

NP 

dt p £ (l n /A n ) 

n 

3 

where q is the volumetric flow rate (m /s) of the propellant, f is a 
coefficient chosen to approximate the Siam of the turbulent line 
losses and the nozzle losses in the feed system (m , 1 and A are 
the length (m) and area (m 2 ) of pipe n, respectively, is the 
pressure in the propellant storage tank, P £ is the pressure in the 
combustion chamber and NP is the total number of pipes in the feed 
system. 

Injector calculations included determination of the droplet 
size distribution and velocity upon injection. Relying heavily upon 


empirical correlations, Webber found the droplet size distribution 
and velocities associated with a self -impinging doublet injector. 

In the combustion chamber, the liquid and gaseous states were 
treated separately. Liquid propellants were divided into droplet 
size groups whose diameters and velocities were followed over time. 
Vaporization rate per droplet, (o v (kg/sec per droplet), was given as 
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The gas state in the combustion chamber was considered well 
mixed and to be a perfect gas. Gas velocity was allowed to vary in 
the combustion chamber. Gas behavior was determined by conservation 
of mass in the combustion chamber, the perfect gas law and an 
expression for the expansion of gas through a choked nozzle. 

Webber integrated his equations numerically using a fixed time 
step integration. The size of the time step was seen to be 
relatively unimportant. He was able to predict a 66 Hz chug whose 
maxi mum amplitude was about 200 psia for an average chamber pressure 
of 400 psia. From his study, Webber concluded that: 

1. chug is sensitive to predicted temperature (he used 
empirical data to give chamber temperature as a function of 
stoichiometry) and 

2. small variations in the droplet size distribution were not 
very important, but gross changes drastically changed the 
model's predictions. 

Calculated frequencies were within 7% of experimental values 
and amplitudes were within 57 .. 

Courtney (1960) used a much simpler approach to determine 
combustion intensity (mass of propellant burned per volume per time 
in a combustion chamber) in a heterogeneous STR. He began with an 
assumed droplet size distribution and assumed steady state 
conditions and infinitely fast chemical kinetics. Droplet 
evaporation rate was given as: 


where k was assumed constant (but not explicitly given), r is the 
droplet radius and <p is the droplet age in the combustion chamber. 
Overall evaporation rate was determined by integrating [2.15] for 0 
< r < r maY . For a monodisperse spray, Courtney showed that 
combustion intensity is a strong function of droplet size; intensity 
increases rapidly with decreasing droplet size. 

2.7 SUMMARY 

Low frequency combustion instability in a liquid bipropellant 
rocket engine involves a coupling of oscillations in chamber 
pressure with oscillations in propellant flow rate. This coupling is 
possible because of the combustion time lag in the combustion 
chamber and the inertia of the propellant in the feed lines. 

A number of models have predicted the domain of operating 
conditions in which a rocket engine will undergo stable operation. 
(Summerf ield, Crocco, Lim and others). These models involve 
linearized equations from which either stability or instability of 
an engine can be predicted. In general, it has been shown from these 
models that higher propellant temperatures, lower chamber pressures 
and longer propellant feed lines lead to more stable rocket engines. 

Other models have predicted either the steady state operation 
or transient behavior of rocket engines (Webber and Courtney). These 


types of models generally involve nonlinear equations which are 
solved over time numerically for variables such as chamber pressure 
and temperature, droplet size distribution and combustion rate. The 
advantage of these models is that amplitudes and frequencies of 
oscillations in pressure, temperature and oxidizer mass flow rate 
can be calculated. These models have shown the STR model to be 
sufficient for modelling-rocket engine combustion chambers 
undergoing chug. 

The author has modelled the SSME fuelside preburner as a 
perfectly mixed reactor and has predicted its behavior during engine 
shutdown. Derivations of the governing equations for this model are 
found in Chapter III. 


CHAPTER III 


MODELLING THE SSME FUEL PREBURNER 

The SSME fuel preburner is modelled as a multi-pipe, multi-node 
oxidizer feed system, a transient stirred tank reactor (TSTR) 
combustion chamber and an exit turbine, as shown in Figure 3.1. This 
approach is similar to that used by Webber (1972) and employs most 
of the same governing equations as those used by Pratt and 
Radhakrishnan (1986). Using a TSTR model allows the frequency and 
amplitude of oscillations in chamber pressure and oxidizer mass flow 
rate to be calculated. Using a multi-pipe, multi-node oxidizer feed 
system allows the position of a liquid/ vapor interface in the feed 
system to be calculated at each integration time step. 

The governing equations for the TSTR and oxidizer feed system 
are partially derived below and their solution is described. Full 
derivations are given in Appendix A. 

The governing equations were solved on a digital computer. A 
program flowchart, listing and sample input and output files are 
found in Appendix B. 

3.1 THE COMBUSTION CHAMBER 

The TSTR model for the combustion chamber assumes that at any 
given time the combustion chamber can be considered a stirred tank 
reactor (properties and composition are uniform throughout the 
combustion chamber). With time, though, chamber composition and 
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Helium Storage Tank (Constant Pressure and Temperature) 



Figure 3.1: Schematic of the Fuelside Preburner Model 



properties are allowed to vary. Using a TSTR to model a chugging 
rocket combustion chamber is justified, in that oscillations in the 
combustion chamber have long periods in comparison to the amount of 
time it takes a disturbance to travel through the combustion 
chamber. Also, it is aggregate, not local changes in pressure in the 
chamber that~ are able to couple with oscillations in the propellant 
feed rate to produce a low frequency system instability. 

The TSTR model was used because it is simple, requiring a 
minimal number of calculations to determine the conditions at each 
time, because it can accommodate heterogeneous combustion and 
because a coupling of oscillations in feed system propellant flow 
with chamber pressure oscillations is possible. 

Heterogeneous combustion of liquid oxygen with gaseous hydrogen 
occurs within the fuel preburner combustion chamber. This is 
accommodated by the TSTR model by the imposition of a uniform 
oxidizer droplet distribution throughout the combustion chamber. 
Droplets were assumed to all have the same diameter and an equal 
number of droplets are found in any part of the combustion chamber 
volume. The instantaneous average droplet diameter (D^) was 
calculated using a relation given by Hersch and Rice (1967). For a 
coaxial injector with gaseous hydrogen injected through an anulus 
around a central stream of liquid oxygen. 


D = c 


m 


dm 




^ 2^2 


Pj T ] 


[3.1] 


H, 
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c, is a dimensionless constant equal to 0.485. (0/F) is the mass 
dm 

ratio of oxidizer flow to fuel flow, is the molecular weight of 

hydrogen (kg/kgmole) , is the area of the anulus through which 

hydrogen is injected (m 2 ), R is the universal gas constant 

8 

(J/kgmole*K), p^ is the density of the oxygen jet (kg/m ), T H is 
the temperature of the hydrogen on injection (K) and P c is the 
chamber pressure (Pa). Equation [3.1] gives the mean droplet 
diameter of liquid oxygen droplets immediately after injection. In 
reality, the mean droplet diameter in the combustion chamber is less 
than that just after injection, however, equation [3.1] was used to 
give the mean droplet diameter throughout the chamber since no data 
were available for determination of the droplet distribution in the 
chamber. The effect that mean droplet diameter has on prebumer 
operation is shown in Chapter IV. 

The evaporation rate of a single droplet, cj i (kg/sec/droplet) , 
was given by Webber (1972) in equation [2.12]. For an average liquid 
density of 3 (mass of liquid per unit volume of the combustion 
chamber), the total vaporization rate, Q (kg/s), including the 
contributions from all droplets, is 


Q = 


6 3V 


P u D 3 
' 1 m 


cu. 


[3.2] 


V is the combustion chamber volume (m ). 


Combustion rate was calculated in two different ways. First, 
all the oxygen that vaporized was assumed to combust immediately 


r 


with hydrogen to form water vapor. This implies that diffusion of 

oxygen vapor into the chamber contents and chemical reaction both 

occur instantaneously. Since oxygen is consumed as soon as it 

vaporizes, its rate of production by combustion T\. (kgmoles/s) is 

u 2 

equal to negative its rate of production by vaporization: 


Q 



From conservation of atoms. 


[3.3] 



For the second combustion model, atomization was assumed infinitely 
fast. Oxygen vaporized according to [2.12] then instantaneously 
mixed with the gaseous species in the combustion chamber. The 
production of water vapor from hydrogen and oxygen occurred via a 
chain reaction involving nine elementary reactions and the species 
H 2 , 0 2 , H 2 0, H, 0 and OH. 

To aid in the solution of the governing equations, the very 
reactive species H, 0 and OH were assumed to be consumed as soon as 
they were produced. In other words, their mole numbers, o^ (kgmoles 
of species i per total mass in the combustion chamber) were constant 
at any time: 
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This condition fixed the values of the o^, Oq and o^ at each time 
step. A full derivation of reaction rate involving chemical kinetics 
is found in Appendix A. 

The governing equations for the combustion chamber were derived 
from the conservation of mass for the gaseous and liquid contents of 
the combustion chamber, conservation of species, conservation of 
energy and the ideal gas relation. 

Conservation of mass for the gaseous contents of the chamber is 
shown in Figure 3.2 and described analytically by equation [3.4]. 


inflow 


outflow 


NSTRM ^ 

£ m. 
Jg 

J-l 


d(p c v ) /dt , Q 


Figure 3.2: Conservation of Mass for the Gas Phase 


d(p c V) NSTRM 


= ( E m. ) - m + Q 
Jg g 


[3.4] 


m is mass flow rate of gases out of the combustion chamber (kg/s), 
g 
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yO is combustion chamber density (kg/m ), t is time (s), V is volume 
3 

(m ), NSTFM is the number of inlet streams to the TSTR, the 

a 

subscript g denotes the gaseous state and superscript denotes an 
inlet quantity. 

Conservation of mass for the liquid contents of the chamber is 
shown in Figure 3.3 and described analytically by equation [3.5]. 


inflow 


NSTRM.* 

E m.. 

j-i Jl 



d(PV)/dt, -Q 


outflow 


Figure 3.3: Conservation of Mass for the Liquid Phase 


d(0V) NSTEM.* 3 

= E m . t - Q - m 

at j-i Pc * 


[3.5] 


The subscript 1 denotes the liquid phase. Since the liquid density 
is uniform throughout the combustion chamber, the liquid outflow is 
given by Pm Ifi . 

Conservation of species is shown in Figure 3.4 and described 
analytically by equation [3.6], 
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inflow 


TSTR 


outflow 


nstrm.* * 

E m.o. . 
3 ij 

j=l 




d(o.p c v)/dt, r. 


►. 

m o. 
g i 


Figure 3.4: Conservation of Species 


d(o.pv) NSTRM * 

— — = E ( m.o. .) - m o. + I\, i=l,2,3, . . .NS [3.6] 

3 ij g i i 

j=l 

o^ is the mole number of species i (kgmoles of i per kg in the 
combustion chamber), NS the number of species present in the 
combustion chamber and I\ the rate of production of species i by 
combustion (kgmoles/s). 

Conservation of energy is shown in Figure 3.5 and described 
analytically by equation [3.7]. 


inflow 

a > 

NSTRM NS .* * * 

E E m.o. .h. . 

J iJ ij 

j=l i=l 



TSTR 

dE/dt 


outflow 

► 

NS 

E m o.h. . 
8 i ij 

i=l 


Figure 3.5: Conservation of Energy 




„ NS NSTRM . . * NS 

AT? • H A • • 

— = E £ m.o. .h. . - £ m o.h. - Q [3.7] 

dt . . . . 3 1J 1J ... 8 1 l 

i=l j=l i=l 

E is the total (chemical plus sensible) energy (Joules) in the 
combustion chamber , h^ is the molar enthalpy of species i 
(J/kgmole) and Q is the heat transfer rate out of the combustion 
chamber (Joules/s). 

Liquids in the chamber were assumed to occupy negligible 
volume. The volume occupied by liquid in the combustion chamber 
after a long time at a steady oxidizer mass flow rate of 5 kg/s and 
a chamber temperature of 1000K was calculated. At these conditions, 
the ratio of the volume occupied by liquid to the total chamber 
volume was 1 : 47,000. Gases were assumed to obey the ideal gas 
relation, which took the form 


P 

c 


NS 

-P c ( .Vi> 

i=l 


R T 
g c 


[3.8] 


NS 

E a. gives the average molecular weight of the chamber gases. 
i=l 1 

Equations [3.6], [3.7] and [3.8] were expressed in terms of the 
logarithmic variables In T , In P £ and In o^, i=l ,2,3, . . .NS. This 
prohibited negative values and facilitated integration by altering 
the time constants associated with these variables. In their final 
form, the conservation of mass for the liquid phase, the 
conservation of species, the conservation of energy and the ideal 
gas relation are 


NSTRM # j. . P 

= ( E m... ) - Q - m 

jlv v gv Q 

i = l > c 


[3.9] 


d(ln o^ 


1 NSTRM* * 

E m . ( o . . - o . ) + T. 

Da. . , JV 1J 1 

rc 1 j=l 

i=l,2,3, . . .NS 


[3.10] 


d(ln T ) 
c 


NS NSTRM >ic * * NS 


E E m. o..h..-Em o.h. - D E u. (da. /dt) 
jv ij ij i=1 gv i l /~c =1 1 1 

i=l j=l 


(dD/dt) E o . u . - Q/V / T D S o.(c . - R ) 
'. =1 ii v cAt i=l 1 P 1 S 


[3.11] 


d(ln P ) 
c 


d(ln T ) 1 dD 1 d NS 

— + — £ + — ( £ a . ) 

dt D dt NS dt i=l 1 

( E a. ) 
i=l 1 


[3.12] 


The superscript v denotes a volumetric quantity (1/m ), c ^ is the 
specific heat at constant pressure for species i (J/kgmole*K) and 
is the universal gas constant ( J/kgmole*K) . The combustion chamber 
equations are derived in greater detail in Appendix A. 


3.2 THE OXIDIZER FEED SYSTEM 


Before shutdown begins, the fuel prebumer oxidizer feed system 


is completely filled with liquid oxygen. During shutdown, gaseous 
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helium progressively replaces oxygen in the pipes. To create a feed 
system which was simple, yet responsive to changes in combustion 
chamber pressure and fluid compressibility in the pipes, the 
following assumptions were made. 

1. Liquid oxygen is incompressible. 

2. Helium behaves as a perfect gas. 

3. Temperature is constant throughout the feed system. 

4. Pipes have a characteristic cross sectional area; area is 
taken as constant throughout a given pipe. 

5. There is no accumulation of mass at pipe nodes; the nodes 
have negligible volume 

6. If a pipe contains both the liquid and the gas phases, a 
distinct interface exists between them. 

7. Fluid velocity is constant within a given pipe (i.e., if 
there is an interface in the pipe, the liquid, the vapor 
and the interface are all assumed to be moving at the same 
velocity. 

8. Flow resistance in a pipe is lumped at the pipe exit, so 
the flow resistance for a pipe containing a liquid/vapor 
interface is the same as that for a pipe which contains 
only liquid, until all of the liquid is expelled. 

Time dependent variables in the feed system include the fluid 
velocity in each pipe, the density at each pipe node and the 
position of the liquid/ vapor interface in each pipe. 

Figure 3.6 shows a pipe containing a liquid/vapor interface. 




Figure 3.6: Schematic of a Pipe Containing a 
Liquid/Vapor Interface 

3 

V is volume (m ), L is pipe length (m), and are upstream and 
downstream velocities (m/s), respectively, and P^ are upstream 
and downstream pressures (Pa), respectively, A is pipe Area (m 2 ) and 
Fj accounts for turbulent and shear head losses in the pipe. 
Subscripts 1 and g denote the liquid and gas phases, respectively. 
Since pipe length is constant, 

dL dL 1 

— & = - — - = v [3.13] 

dt dt 

Conservation of momentum, expressed in integral form, is 
(White, 1979) 
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[3.14] 


/ 


# 


S F =^~(||| P vdV) + JJ P v( v*n) dA 
CV CS 


2 F is the sum of the longitudinal forces on the pipe. From Figure 
3.6, one can see that 


n ■ (Pj ■ P 2 ) A - F f _ [3.15] 


Fj is given as 


F f = f L p v 2 A 


[3.16] 


where f is a constant (1/m) corresponding to the longitudinal force 
that either the gas or the liquid within the pipe is exerting on the 
pipe. Recall that if any liquid is in the pipe, it is as though 
there is only liquid in the pipe. 

Using equations [3.13], [3.15] and [3.16], equation [3.14] can 
be written 


8 dt J 


dp. D L p.L. 

-LI + { J±JL + — = 


dt 


dv 

dt 


P 2 - fL Pv 2 
v 


[3.17] 


Equation [3.17] is applied to each pipe, leading to NP equations (NP 
is the number of pipes in the feed system) . 

Figure 3.7 shows a pipe node. 
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Figure 3.7: Schematic of a Pipe Node. 


Since mass is not allowed to accumulate at the node, for NPN pipes 
connected at the node, conservation of mass for the node can be 
written 


NPN 

1 ( Pn A n V n> = 0 [3.18] 

, < n n n 

n=l 

Equivalently, 

\ 

Taking the derivative in [3.19] and assuming that the fluid in each 
pipe has the same density as the fluid in the others, equal to the 
density in the node, leads to 


NPN . 

Z — GO A V ) = 0 
. . r'n n n 

n=l dt 


[3.19] 


[3.20] 


d A 

dt 


NPN 

2 A v 
, n n 
n=l 


NPN 

n 2 A (dv /dt) = 0 
r'n . n n 
n=l 


[3.20] is applied to NN nodes, where NN is the number of nodes in 
the oxidizer feed system. The piping system equations are derived in 
greater detail in Appendix A. 

Equations [3.17] and [3.20] constitute NN+NP equations in the 
NN unknown nodal densities and the NP unknown pipe velocities. This 
set of equations is linear and, at each time step, can be solved for 
dv^/dt (n=l,2,3,... NP) and d/^/dt (k=l,2,3,... NN) . 

3.3 THE EXIT TURBINE 


Hot product gases pass through a turbine when exiting the fuel 
prebumer. The mass flow rate of the gases through the turbine is a 
function of the pressure downstream of the turbine (PHG) and the 
inertia of the turbine blades, as well as the pressure in the 
combustion chamber. An equation which gives the mass flow rate 
through the turbine was given in the Dynamic Balance Model for the 
SSME (Nguyen, 1981) as 
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The parameter (c. ) is an empirical function of turbine geometry 

ftp 

and speed. For this model, c^ i s assumed to be constant during the 
short time duration of the chug. This implies that turbine speed 
decreases slowly. 


3.4 SOLUTION OF THE GOVERNING EQUATIONS 

Equations [3.9], [3.10], [3.11], [3.12], [3.17] and [3.20] were 

integrated numerically to give the values of T , P , 

(i=l,2,3,...NS),p (k=l,2,3, . . .NN), v. ( j=l,2,3, . . .NP) and L 

J S 

(position of the liquid/vapor interface in pipe j) ( j=l,2,3, . . .NP) 
over time. The TSTR equations are nonlinear in the reaction rate 
expressions and the expression that give thermodynamic properties. 
Since the set of equations is also stiff (involving widely varying 
time constants), the IMSL math subroutine DGEAR was used to 
integrate the equations. DGEAR finds approximations to the solution 
of ordinary differential equations of the form 

y’ = f(x,y) 

for given initial conditions. The function f(x,y) was calculated in 
a user supplied subroutine which was declared external in the main 
program and was accessed by DGEAR for each time step in the 
integration. A number of different integration methods were 
available in the DGEAR package. For most executions, the stiffness 
methods of Gear were used. Where possible, initial conditions were 
chosen to agree with experimental data. In other cases, when 


experimental data were not available, "steady state" conditions were 
used as initial conditions. The selection of initial conditions is 
further discussed in Chapter IV. 

Flow charts showing the solution process, a listing of program 
TRNCHG which was developed to solve the governing equations using 
DGEAR, a sample input file and sample output files are found in 
Appendix B. 
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CHAPTER IV 


RESULTS 

In Chapter III a Transient Stirred Tank Reactor (TSTR) model 
was proposed for predicting the transient behavior of the SSME 
fuelside preburner during the helium purge of the engine shutdown. 
Fuel feed rate is assumed constant, while the oxidizer feed system 
is modelled as a series of pipes and nodes and an injector 
resistance. The exiting fluid from the combustion chamber, 
containing gaseous products, gaseous reactants and liquid droplets, 
passes through an exit turbine. 

A computer program, TRNCHG (Appendix B), was written to 
numerically integrate the resulting differential equations. From the 
results of TRNCHG it was possible to calculate the amplitude and 
frequency of oscillations in combustion chamber pressure and 
oxidizer mass flow rate. A number of engine parameters were varied 
to determine the role they play in the fuelside preburner chug and 
to demonstrate the reliability of TRNCHG. The parameters that were 
varied include the rate of heat transfer from the walls of the 
combustion chamber to the contents of the chamber, the mean droplet 
diameter, the fuel mass flow rate, the fuel temperature, the liquid 
oxygen (LOX) temperature, the temperature of the helium, and the 
length of the helium purge line (the pipe connecting the helium 
storage tank to the helium check valve (Figure 3.1, p. 43)) and the 
diameter of the helium purge line. 


This chapter is a summary of the results of TRNCHG. All marked 
data points in the figures that follow were generated mathematically 
and were marked either to differentiate between curves on the same 
graph or to specify calculation points on curves drawn using few 
data points. 

4.1 OPERATION OF THE FUELSIDE PREBURNER DURING THE SSME SHUTDOWN 
HELIUM PURGE 

As described in Chapter I, during the SSME shutdown, gaseous 
helium pushes liquid oxygen from the oxidizer feed system of the 
fuelside preburner and into the combustion chamber. As outlined in 
Table 1.1 (page 7), just before the helium purge of the fuelside 
preburner begins, the fuelside preburner oxidizer valve (FPOV) is 
closed, so the velocity of the LOX in the prebumer oxidizer feed 
system was set at 0.001 m/s (it was not set equal to zero because of 
the numerical difficulties arising from setting the velocity equal 
to zero). It is possible that a small amount of LOX enters the 
preburner combustion chamber between the time that the FPOV is 
closed and the valve that restrains the helium is opened. This is 
due to the heating and expansion of LOX in the pipe connecting the 
FPOV to the preburner manifold. 

The pressure throughout the combustion chamber, oxidizer feed 
system and helium supply piping was assumed uniform and equal to the 
cracking pressure of the helium check valve, 5.17 MPa (750 psia). 


i 


The NASA SSME transient model predicts that the temperature in 
the fuelside prebumer decreases at a rate of about 50 K/s (90°R/s) 
(Seymour, 1986) at the beginning of shutdown. Without any heat 
transfer from the chamber walls, TRNCHG predicts that the fuelside 
preburner temperature falls at a rate of 280,000 K/s when the helium 
check valves first open. Because the rate at which TRNCHG predicts 
the chamber temperature decays was too high, a rate of heat transfer 
to the chamber contents of 6.5 x 10^ J/s (6.2 x 10^ Btu/s) was 
selected for the TRNCHG helium purge simulation. This value was 
chosen because it made the initial temperature gradient predicted by 
TRNCHG closer to that predicted by the SSME transient model yet did 
not fix the chamber temperature. The selection of heat transfer rate 
will be further discussed in section 4.2. 

At the beginning of the fuelside preburner helium purge, the 
pressure downstream of the exit turbine (PHG) is 2.41 MPa (350 
psia). By the time the helium purge is over, PHG is 1.03 MPa (150 
psia). The average value of PHG, 1.72 MPa (250 psia), was used for 
the helium purge simulation. The exit mass flow rate coefficient, 
c^ , was calculated using equation [3.21] (page 56). c^ was 
algebraically isolated and calculated by setting the initial exit 
mass flow rate equal to the fuel mass flow rate (which is constant 
during shutdown) and setting the chamber pressure equal to the 


initial chamber pressure, 5.17 MPa (750 psia). 


c-. was assumed 
ftp 


constant during the short time duration of the chug. 
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No data on helium temperature in the helium line were 
available, so the temperature of the helium purge was set equal to 
the oxidizer temperature, 120 K (216°H). 

A summary of inputs to TRNCHG for the SSME shutdown is given in 
Table 4.1. 

Figure 4.1 shows TRNCHG' s prediction for combustion chamber 

temperature during the helium purge. Time equal to zero corresponds 

to the time that the helium check valves are opened for this graph 

and all graphs generated for the SSME shutdown simulation. Chamber 

temperature falls until it reaches a fairly constant value, around 

400 K (720°R). The steady value that the temperature reached was 

determined largely by the rate of heat transfer to the chamber 

contents. With no heat transfer, the steady state temperature is 188 

K (338°R). Figure 4.1 includes only the first 25 msec of TRNCHG' s 

predictions although the program was run until the helium/LOX 

interface in the oxidizer feed system reached the combustion 

chamber. Figure 4.2 is a plot of the mole number of oxygen vapor 

(kgmoles of oxygen vapor per mass (kg) in the combustion chamber) in 

the combustion chamber during shutdown versus time. The amount of 

oxygen vapor in the chamber rises sharply (from the program's 

-12 

minimum allowable value of 1.0 x 10 kgmoles/kg) to a steady value 

_3 

around 5 x 10 kgmoles/kg. The mole number of water vapor did not 
change significantly from the minimum allowable mole number. 


Table 4.1: Inputs to TRNCHG for the SSME Shutdown Simulation 


PARAMETER 

VALUE AT THE BEGINNING 
OF THE SSME 
SHUTDOWN SIMULATION 

Chamber Temperature 

550 

K 

(990°R) 

Chamber Pressure 

5.17 

MPa 

(750 psia) 

PHG 

1.72 

MPa 

(250 psia) 

LOX Temperature 

120 

K 

(216°R) 

Fuel Temperature 

160 

K 

(288°R) 

Helium Temperature 

120 

K 

(216°R) 

Fuel Mass Flow Rate 

21.0 

kg/s 

(46.3 lbm/s) 

Heat Transfer Rate 

6.5 x 10 7 

J/s 

(6.2 x 10^ Btu/s) 

Pipe Length from the 

1.5 

m 

(4.9 ft) 

Helium Storage Tank 




to the Check Valve 
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Figure 4.1: Temperature v. 


Time for SSME Shutdown 
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Figure 4.2: Mole Number of Oxygen Vapor in the Combustion 
Chamber v. Time for SSME Shutdown 
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As shown by Figures 4.1 and 4.2, TRNCHG predicts that ignition 
will not take place in the SSME fuelside preburner during the 
shutdown helium purge, contrary to the findings of test firings. 
There are two explanations for TRNCHG's incorrect prediction. First, 
using a TSTR model results in the overprediction of the equivalence 
ratio at which combustion takes place. In the heterogeneous stirred 
tank reactor, LOX vaporizes and the oxygen vapor is assumed to 
instantaneously mix with the gaseous contents of the combustion 
chamber. In actual droplet combustion, though, oxidizer and fuel 
burn on a diffusion flame surface around the droplet at a position 
where oxidizer and fuel meet in stoichiometric proportion. Thus, the 
TSTR model overpredicts the equivalence ratio at which combustion 
takes place; the mixture of gases is beyond the rich limits of 
combustion. Second, the kinetic mechanism and elementary reaction 
rate data used in TRNCHG may be inaccurate at the temperatures and 
pressures that are found in the preburner during the helium purge. 

Since use of chemical kinetics in the calculation of combustion 
rate resulted in erroneous predictions, a second combustion model 
was used to predict preburner behavior during shutdown. All the LOX 
that was vaporized was assumed to immediately combust with hydrogen 
to form water vapor (i.e., chemical kinetics were assumed to be 
infinitely fast). Figure 4.3 shows TRNCHG's prediction of chamber 
pressure using this second combustion model and Figure 4.4 shows the 
predicted oxidizer mass flow rate. The symbols mark each third 
analytically generated data point. TRNCHG was allowed to run until 
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Figure 4.3: Chamber Pressure v. Time for the SSME Shutdown 
(Infinitely Fast Chemical Kinetics) 
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Figure 4.4: Oxidizer Mass Flow Rate v. Time for the SSME Shutdown 
(Infinitely Fast Chemical Kinetics) 
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the helium/LOX interface entered the combustion chamber, but only 

the first 25 msec of the run is shown in Figures 4.3 and 4.4, since 

no oscillations took place after that time. The initial peak to peak 

amplitude of the pressure oscillation is 1.4 MPa (200 psia) and the 

frequency of the oscillation is 111 Hz. The amplitude and frequency 

of pressure oscillations during an actual typical preburner chug are 

2.41 MPa (350 psia) and 117 Hz, respectively. Oscillations in the 

chamber pressure and oxidizer mass flow rate are about 1/5 cycle 

more than 180° out of phase, due to the vaporization time delay. 

To investigate the possibility that a disturbance in chamber 

pressure (for example, due to a change in the oxidizer preburner 

combustion chamber pressure) incites the fuelside preburner chug 

when the helium/LOX interface has progressed nearly to the 

combustion chamber, a perturbation was applied in TRNCHG at 50 msec. 

The system was perturbed by changing the constant, c^, in the 

expression for mean droplet diameter (equation [3.1], page 44). c^ 

was changed from 0.485 to 20.0 at time equal to 50 msec. After 

2.5 msec the value of c, was returned to 0.485. While c, was equal 

dm dm 

to 20.0, the vaporization rate decreased since the droplets were 

larger and the ratio of the surface area of liquid to the mass of 

liquid in the combustion chamber decreased. Because of the 

subsequent decrease in the evaporation rate of the LOX, liquid 

accumulated in the combustion chamber. When the value of c, was 

dm 

returned to 0.485, the accumulated liquid burned quickly and the 
chamber pressure rose sharply. The resulting chamber pressure 


oscillations are shown in Figure 4.5. Once again, chamber pressure 
oscillations damped quickly. 

Assuming infinitely fast chemical kinetics, TRNCHG did not 
predict sustained oscillations in chamber pressure. As described 
earlier, a number of the engine parameters used in this base case 
simulation of the SSME shutdown were selected as best guesses of the 
actual conditions that exist in the fuelside preburner during 
shutdown. To determine whether or not a change in engine parameters 
could incite a chug, engine parameters were varied, one at a time. 

In many cases a small change in an engine parameter had a large 
influence on preburner stability. Following is a summary of the 
predictions of TRNCHG when engine parameters were varied. 

4.2 THE EFFECT OF HEAT TRANSFER RATE ON PREBURNER 
OPERATION 

As discussed in section 4.1, in order for TRNCHG' s prediction 
of initial temperature gradient to be equal to that predicted by the 
NASA SSME transient model, there must be heat transfer to the 
combustion chamber contents. Care had to be taken in choosing this 
rate of heat transfer, since data from which heat transfer rate 
could be calculated was not available and since heat transfer rate 
has such a large effect on the chamber temperature. 

TRNCHG was executed using different rates of heat transfer and 
no oxidizer mass flow to find the temperature that the chamber 
contents would eventually reach. These steady state temperatures are 
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Figure 4.5: Prebumer Response to a Perturbation During Shutdown 
(Infinitely Fast Chemical Kinetics) 



Mean droplet diameter proved very important in determining the 

operation of the fuelside preburner, to vary the mean droplet 

diameter, the constant, c^, in equation [4.1], was given the values 

0.0485, 0.97, 1.94 and 3.88. These droplet coefficient values 

resulted in mean droplet diameters (averaged over the 25 msec shown 

in Figure 4.8) of 1.2 ym, 24.4ym, 63ym and 440ym, respectively. Note 

that changing c^ m changed the droplet diameter both directly and 

indirectly, since different values of c^ lead to different values 

of average chamber pressure and oxidizer mass flow rate. Figure 4.8 

shows predicted chamber pressure for the four values of the droplet 

coefficient and Figure 4.9 shows oxidizer mass flow rate. For c, 

am 

below 1.0, pressure oscillations damp quickly. For c^ greater than 
1.0, oscillations become more sustained, until droplet size becomes 
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Figure 4.6: Steady State Temperature v. Heat Transfer Rate 
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CHAMBER PRESSURE (MPa) 



Figure A. 7: The Effect of Heat Transfer Rate on Chamber 
Pressure Oscillations 
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Figure 4.8: The Effect of Mean Droplet Diameter on Chamber 

Pressure Oscillations 
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so large that vaporization rate becomes very small (due to limited 

droplet surface area) and reaction blows off, as shown by the curve 

for c^ m = 3.88. For c^ = 0.0485 the frequency of the pressure 

oscillation is 116 Hz. For c^ = 1.97 the frequency is 83 Hz. For 

c, =1.97 backflow of combustion chamber contents into the oxidizer 
dm 

feed system is predicted, as shown in Figure 4.9. 

As seen in Figure 4.8, there is probably some critical mean 

droplet diameter (or equivalently, some critical combustion delay 

time) at which the fuelside preburner begins to undergo unstable 

operation. The critical diameter could be reached by changing the 

injector geometry or the fuel mass flow rate. The effect of fuel 

mass flow rate on the fuelside preburner is described in section 

4.5. Since, in the absence of chemical kinetics, the value of c, 

dm 

determines the combustion time delay, c^ could be raised to 

increase the combustion time delay and account for chemical 

kinetics. The exact amount that c, would have to be raised was not 

dm 

calculated by the author. 

4.4 THE EFFECT OF THE OXIDIZER FEED SYSTEM ON PREBURNER OPERATION 

TRNCHG was executed for LOX temperature equal to 90 K, 100 K, 
120 K and 150 K (162°R, 180°R, 216°R and 270°R). All the other 
engine parameters were the same as those given in Table 4.1 (page 
64). Decreasing the LOX temperature increases the LOX density and 
decreases the LOX vaporization rate in the combustion chamber. 
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Predicted chamber pressure is plotted against time in Figure 4.10. 
As predicted by Lim (1986), lower oxidizer temperatures resulted in 
less stable preburner operation (i.e., larger amplitude pressure 
oscillations that damped more slowly). 

The frequency of chamber pressure oscillations is plotted 
against LOX temperature iri~~Figure 4.11. As LOX temperature is 
decreased, oscillation frequency decreases slightly. 

In this study, the term amplitude ratio refers to the ratio of 
the peak to peak amplitude of the second chamber pressure 
oscillation to the peak to peak amplitude of the first oscillation. 
It is used as a measure of the rate at which pressure oscillations 
damp out. Amplitude ratio is plotted against LOX temperature in 
Figure 4.12. The amplitude ratio rises linearly as LOX temperature 
is decreased from 150 K to 100 K. As LOX temperature is decreased 
below 100 K, the amplitude ratio increases more slowly. 

4.5 THE EFFECT OF THE FUEL FEED SYSTEM ON PREBURNER OPERATION 

Two fuel feed system parameters, the fuel temperature and the 
fuel mass flow rate, were varied in TRNCHG. As predicted by Lim 
(1986), lower fuel temperature led to sustained pressure 
oscillations. Contrary to Lim's predictions, fuel mass flow rate 
also influenced the preburner's stability. 

TRNCHG was executed for fuel temperatures of 50 K, 100 K and 
200 K (90°R, 180°R and 360°R). The resulting chamber pressure was 
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The Effect of Oxidizer Temperature on Chamber 
Pressure Oscillations 
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Figure 4.12: Amplitude Ratio v. LOX Temperature 
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plotted against time in Figure 4.13. Lower fuel temperatures 
resulted in more slowly damped pressure oscillations of greater 
amplitude and lower frequency. Oscillation frequency is plotted 
against fuel temperature in Figure 4.14. Frequency does not change 
much with temperature until the temperature becomes less than about 
120 K, where frequency decreases rapidly as fuel temperature is 
decreased. 

In Figure 4.15, chamber pressure was plotted against time for 
executions of TRNCHG using fuel mass flow rates of 12 kg/s, 18 kg/ s 
and 25 kg/s. So that the results of these executions were not 
obscured by heat transfer rate to the chamber contents, the heat 
transfer rate was assumed to be 0.0 J/s for these runs. The fuel 
mass flow rate did not effect the pressure oscillation frequency, 
but did influence the rate at which the oscillations damped. The 
average oxidizer mass flow rate is plotted against the fuel mass 
flow rate in Figure 4.16. 

The average oxidizer mass flow rate rises almost linearly with 
fuel mass flow rate. For high mass flow rates (greater combustion 
chamber throughput) , more energy is needed to heat the fuel that 
does not participate in combustion and hot product gases are quickly 
swept from the combustion chamber. This leads to lower chamber 
pressures which lead to higher oxidizer mass flow rates. 

As mentioned above, Lim (1986) predicted that fuel mass flow 
rate is not important in the fuelside prebumer chug. The unexpected 
behavior of chamber pressure oscillations as mass flow rate is 
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Figure 4.13: The Effect of Fuel Temperature on Chamber 
Pressure Oscillations 
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Figure 4.14: Pressure Oscillation Frequency v. Fuel Temperature 
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Figure 4.15: The Effect of Fuel Mass Flow Rate on Chamber 
Pressure Oscillations 
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Figure 4.16: Oxidizer Mass Flow Rate v. Fuel Mass Flow Rate 
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changed is attributed to oscillations in chamber temperature. Recall 
that the fuels ide preburner operates very fuel rich and that the 
incoming fuel is quite cold (160 K [288°R]). For a high fuel mass 
flow rates, much of the energy liberated by combustion raises the 
temperature of the hydrogen (most of which does not participate in 
chemical reactions) to the chamber temperature; there is not enough 
energy to raise the temperature much above the chamber temperature 
and sustain oscillations. For a lower mass flow rate, less energy is 
needed to heat the incoming hydrogen to the chamber temperature, so 
an oscillation in chamber temperature and an oscillation in chamber 
pressure can be sustained. 

4.6 THE EFFECT OF THE EXIT TURBINE ON PREBURNER OPERATION 

During shutdown, the pressure downstream of the exit turbine, 
PHG, falls from 2.41 MPa (350 psia) to 1.03 MPa (150 psia). It is 
possible that this change plays a part in inciting the fuel 
preburner chug. Figure 4.17 contains plots of chamber pressure 
versus time for PHG equal to 1.03 MPa, 1.72 MPa and 2.41 MPa. 
Frequency of the pressure oscillations is the same for all three 
values of PHG. The rate at which the oscillations are damped is 
smaller for lower values of PHG, as shown in Figure 4.18. Although 
changing the pressure downstream of the exit turbine did not produce 
chugging by itself, this change probably plays a part in the onset 



Figure 4.17: The Effect of the Pressure Downstream of the Exit 
Turbine on Chamber Pressure Oscillations 


♦ 


88 




4.7 THE EFEECT OF THE HELIUM PURGE ON PREBURNER OPERATION 


At the onset of this research, it was thought that oscillations 
in helium density in the oxidizer feed system were at least partly 
responsible for the occurrence of the fuelside preburner chug. In 
section 4.1 it was shown that, even when the oxidizer feed system 
was almost completely filled with helium, a perturbation did not 
result in chug, so the presence of helium does not single-handedly 
bring about chug. To see what effect the presence of helium does 
have on fuel preburner operation, the temperature of the helium 
(assumed to be uniform and constant) and the helium purge line 
diameter and length were varied. 

No data were available for the temperature of the helium at the 
helium check valve, so the temperature was initially assumed to be 
equal to the LOX temperature, 120 K (216°R). In Figure 4.19, chamber 
pressure is plotted against time for helium temperatures of 50 K, 

150 K, 200 K and 300 K. Figure 4.20 shows the variation of the 
frequency of chamber pressure oscillations with helium temperature 
and Figure 4.21 shows the variation with helium temperature of the 
ratio of the peak to peak amplitudes of the second and first 
pressure oscillations. Higher helium temperatures lead to higher 
frequency pressure oscillations that damp more slowly. In the 
absence of helium, higher temperatures in the oxidizer feed system 
lead to more stable operation. Helium temperature reverses this 
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Figure 4.19: The Effect of Helium Temperature on Chamber 
Pressure Oscillations 
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Figure 4.21: Amplitude Ratio v. Helium Temperature 
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trend. The presence of helium in the oxidizer feed system results in 
less stability for higher temperatures in the oxidizer feed system. 
The initial amplitude of pressure oscillations differed only 
slightly for the different helium temperatures. 

As predicted by Summerf ield (1951), shorter pipe lengths lead 
to chug. Figure 4.22 is a plot of chamber pressure versus time for 
the fuels ide pre burner for no helium purge line (i.e., for the 
helium tank mounted directly on the helium check valve). Execution 
of TRNCHG was halted after 8.3 msec because execution became very 
slow. For no helium line, the preburner undergoes undamped 
oscillations in chamber pressure and backflow of chamber contents 
into the oxidizer feed system. 

Figure 4.23 is a plot of chamber pressure versus time for 
helium purge line lengths of 0.5 m and 2.0 m. Oscillation frequency 
and the ratio of the peak to peak amplitude of the second 
oscillation to that of the first are plotted against helium line 
length in Figures 4.24 and 4.25, respectively. Both the frequency 
and amplitude ratio rise sharply for helium line lengths below 
0.5 m. 

TRNCHG was executed using helium purge line diameters of 
1.029 cm, 0.739 cm, 0.381 cm and 0.178 cm (0.4", 0.291", 0.381" and 
0.15"). The resulting plots of chamber pressure versus time are 
found in Figure 4.26. For a helium purge line diameter of 0.178 cm, 
TRNCHG was unable to dependably calculate chamber pressure after 
about 5 msec because of numerical instability. The velocity of the 
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Chamber Pressure v. Time, No Helium Purge Line 
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4.23: The Effect of Helium Line Length on Chamber 
Pressure Oscillations 
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Figure 4.25: Amplitude Ratio v. Helium Line Length 
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Figure 4.26: The Effect of Helium Line Diameter on Chamber 

Pressure Oscillations 


fluid in the purge line began to undergo undamped high frequency 
numerical oscillations. To try to prevent numerical instability, the 
maximum error allowable in the IMSL integration routine, DGEAR, was 
reduced and the number of pipes in TRNCHG that represented the purge 
line was increased. Both of these changes resulted in slightly 
longer periods of numerical stability accompanied by considerably 
longer CPU times. The CPU time necessary for TRNCHG to give an 
accurate solution for a small diameter helium purge line was deemed 
too high and only the limited solution shown in Figure 4.26 was 
generated for a helium purge line diameter equal to 0.178 cm. 

From Figure 4.26, one can see that a larger diameter helium 
purge line results in higher frequency chamber pressure oscillations 
of greater amplitude and which damp more slowly. The average chamber 
pressure decreases with decreases helium purge line diameter. For a 
given helium line length, the greater the volume of helium in the 
line, the less stable the preburner. 

4.8 PREBURNER OPERATION AT FULL POWER LEVEL 

TRNCHG was used to simulate fuelside preburner operation at 
full power level (with no helium purge) both to validate the program 
and to determine the part that chemical kinetics play in the 
prebumer chug. Care had to be taken in choosing TRNCHG' s initial 
conditions, since poor choices resulted in large derivatives and 
numeric overflows. Species mole numbers and the initial mass of 


t 


liquid in the combustion chamber were set equal to the values they 
would reach in a transient stirred tank reactor, after a long period 
of time, at constant chamber pressure, chamber temperature and 
oxidizer and fuel mass flow rates. The oxidizer mass flow rate was 
set equal to the oxidizer mass flow rate that would be reached in 
the pipes afte~^ a long period of time at constant upstream and 
chamber pressures. A summary of inputs to TRNCHG for the full power 
level simulation is given in Table 4.2 

TRNCHG was first run assuming infinitely fast chemical kinetics 
and was then run taking chemical kinetics into account (as described 
in Chapter III and Appendix A) . Chamber pressure for these runs is 
plotted against time in Figure 4.27. Just as in actual operation, 
the preburner does not undergo sustained chug at full power level. 
The initial pressure oscillations are due to the choice of the 
initial values of chamber pressure and chamber temperature. The 
presence of chemical kinetics has a pronounced effect on chamber 
pressure and temperature: pressure oscillations damp more slowly and 
the temperature and pressure that the combustion chamber reaches 
after a long time are lower when chemical kinetics are considered. 
The frequency of the pressure oscillations for the run involving 
kinetics is 164 Hz. 

A perturbation was applied to the preburner, at full power 
level, with chemical kinetics included, to see if chug could be 
incited by a large disturbance. At 30 msec, the mean droplet 
diameter coefficient, c^ m> of equation [4.1], was reduced from 0.485 
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Table 4.2: Inputs to TRNCHG for Full Power Level 


PARAMETER 

VALUE AT THE BEGINNING OF THE 
FULL POWER LEVEL SIMULATION 

Chamber Temperature 

1106 K (1991°R) 

Chamber Pressure 

38.0 MPa (5500 psia) 

Mass of Liquid in the Chamber 

2. 71xl0~ 2 kg (5.01xl0~ 2 lbm) 

Oxidizer Mass Flow Rate 

32.08 kg/s (70.72 lbm/s) 

Fuel Mass Flow Rate 

34.40 kg/s (75.84 lbm/s) 

Pipe Length 

0.12 m (0.39 ft) 

Mass Fractions: 


He 

0.0 

H 2 

5.40xl0 _1 

°2 

6.59xl0' 3 

h 2 o 

4.54X10' 1 

H 

5.76xl0' 6 

0 

1.53xl0' 7 

OH 

1.15xl0" 7 






N 

fO 


I • ' 

0 


• i ' 

10 


i — i — | i — i i i — | — i i I i | i r 

20 30 40 

TIME (msec) 


Figure A. 27: Chamber Pressure v. 
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to 0.364 for 0.5 msec, and then returned to 0.485. The resulting 
chamber pressure oscillation is shown in Figure 4.28. Initially, 
when c ^ m is reduced, the mean droplet size becomes smaller and more 
oxygen vaporizes and is available for combustion. During this time 
the pressure rises. Next, when the amount of liquid in the chamber 
and the injection rate of LOX into the chamber become small, the 
combustion rate and pressure fall. The oxidizer mass flow rate then 
increases and the oscillation in chamber pressure continues. Within 
seven cycles, the peak to peak amplitude of the pressure 
oscillations is damped to 16% of its original value. 

4.9 CONDITIONS THAT RESULT IN INSTABILITY 

Aside from reducing the helium line length to zero, it was not 
possible to incite chug by varying just one engine parameter, so it 
is a combination of engine parameters that cause the fuelside 
preburner chug. TRNCHG was run with various combinations of engine 
parameters. The purpose of these runs was not to determine the 
domain of engine parameters that give rise to chugging, but to 
demonstrate that TRNCHG can predict chug and to show that only small 
deviations from the engine parameters used for the SSME shutdown 
simulation can result in chug. Two runs that resulted in instability 
are described below. 

Figure 4.29 is a plot of chamber pressure v. time for hydrogen 
temperature equal to 120K, helium temperature equal to 200 K, PHG 
equal to 1.03 MPa, fuel mass flow rate equal to 18 kg/s and mean 


102 


CHAM 




o ~ 
• «*_ 

h- 

n ■ 


o i 

(£> |~i r n~"| rn i | i i i i | i i i i j m r~f~i r ~ r~T~] t i t i | r i i t | » > 

^ 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 

TIME (msec) 


Figure 4.28: Preburner Response to a Perturbation at Full 

Power Level 
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Figure A. 29: Unstable Preburner Operation (a) 



droplet diameter coefficient equal to 0.585. All other engine 
parameters were the same as those listed in Table 4.1. Execution was 
halted after 49 msec because time for each integration step became 
too long. Using a mean droplet diameter of Q.485 and all other 
engine parameters the same as those described for Figure 4.29, 

TRNCHG predicted stable operation. 

Figure 4.30 is a plot of chamber pressure versus time for 
hydrogen temperature equal to 70 K, fuel mass flow rate equal to 15 
kg/s, heat transfer rate to the chamber contents equal to 8.5 x 10^ 
J/s, PHG equal to 1.03 MPa and helium temperature equal to 200 K. 

All other engine parameters were the same as those listed in Table 
4.1. This run was be made stable by lowering the helium temperature 
to 150 K. 

4.10 SUMMARY OF RESULTS 

The program TRNCHG, which was written to solve the governing 
equations derived in Chapter III for heterogeneous preburner 
operation, was supplied with the approximate conditions at the 
opening of the helium check valves and was executed. When chemical 
kinetics were allowed to proceed according to a chemical kinetic 
mechanism, the model erroneously predicted that ignition would not 
take place in the prebumer. When chemical kinetics were assumed 
infinitely fast, the model predicted that a quickly damped pressure 
oscillation with a frequency of 111 Hz would take place immediately 
after the helium check valves were opened. An oscillation in 



Figure 4.30: Unstable Preburner Operation (b) 




oxidizer mass flow rate into the preburner combustion chamber was 
also predicted and it was about 1/5 cycle more than 180° out of 
phase with the pressure oscillation. 

To determine their effect on engine stability, the following 
engine parameters were varied: heat transfer rate from the chamber 
walls, mean droplet diameter, oxidizer temperature, fuel 
temperature, fuel mass flow rate, pressure downstream of the fuel 
preburner exit turbine (fuel turbo-pump), helium temperature, helium 
line length and helium line diameter. It was found that small 
variations in a number of these parameters had a large effect on 
preburner operation. A summary of the effect that each of these 
parameters has on preburner operation is given in Table A. 3. 

TRNCHG was also used to predict conditions in the fuel 
preburner at full power level. This was both to validate the model 
and to allow comparison of different reaction rate models. As 
expected, TRNCHG predicted stable operation at full power level. 
Assuming that rates of production of the species in the combustion 
chamber could be calculated using elementary reaction rates resulted 
in much more slowly damped oscillations in chamber pressure than 
asstuning chemical kinetics were infinitely fast. 

Finally, it was demonstrated that, after changing a small 
number of engine parameters by only a small amount, TRNCHG was able 
to predict unstable preburner operation. 
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Table 4.3: Summary of the Effects of Engine Parameters on 
Fuels ide Preburner Operation 


ENGINE PARAMETER 
BASE VALUE 

RANGE OF VALUES USED IN 
TRNCHG 

EFFECT ON PREBURNER OPERATION 

Heat Transfer Rate 
7 

Heat transfer rate had a large 
effect on average chamber pressure 

Q = 6.5x10 J/s 

0.0 J/s < Q < l.OxlO 8 J/s 

and temperature, but very little 
effect on the rate at which pres- 

sure oscillations damped. 

Mean Droplet Diameter 

Higher values of c^ lead to large 
amplitude pressure oscillations 
which damp more slowly and have 

Coefficient 

c , « 0.485 

dm 

lower frequencies. For some value 
of c^ between 1.94 and 3.88, 
there is a critical value of c^ 
for which the preburner is un- 
stable. Above that value liquid 
droplets become too large and 
ignition does not occur. 

0,0485 < c . < 3.88 

dm 

LOX Temperature 

Lower oxidizer temperatures lead 
to higher frequency chamber pres- 

T = 120 K 

sure oscillations which damp rela- 

OX 

tively quickly. The amplitude 

90 K < T < 150 K 

ratio for T = 150 K is 54.4% of 

ox 

that for T ° X = 90 K. 
ox 

Fuel Temperature 

Below 120 K fuel temperature has 
a large effect on preburner sta- 

T- . = 160 K 
fuel 

bility. The frequency of chamber 
pressure oscillations for 

50 K < T, . < 200 K 
fuel 

T* = 50 K is 68% of that for 

Ifuel = 150 K. Oscillations damp 
more slowly and have greater am- 
plitudes for lower fuel 
temperatures 
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Table 4.3: (continued) 


ENGINE PARAMETER 
BASE VALUE 

RANGE OF VALUES USED IN 
TRNCHG 

EFFECT ON PREBURNER OPERATION 

Fuel Mass Flow Rate 
= 21.0 kg/s 
12.0 < m f < 25.0 

Decreasing the fuel mass flow 
rate caused chamber pressure osci- 
lations to damp more slowly, but 
did not change the oscillation 
frequency. Higher fuel mass flow 
rates resulted in higher oxidizer 
mass flow rates. 

Pressure Downstream of 

Decreasing the pressure down- 

the Exit Turbine 

stream of the exit turbine result- 

PHG =1.72 MPa 

1.03 MPa < PHG < 2.41 MPa 

ed in more slowly damped chamber 
pressure oscillations, but did 
not change the frequency. The am- 
plitude ratio for PHG = 1.03 MPa 
is 4.9 times that for 
PHG = 2.41 MPA 

Helium Temperature 

T„ = 120 K 
He 

50 K < T„ < 300 K 
He 

Changing the helium temperature 
changed both the amplitude ratio 
and the oscillation frequency; 
Higher helium temperatures lead 
to higher frequency pressure 
oscillations that damped more 
slowly. The frequency at 200 K is 
1.3 times that at 120 K. The am- 
plitude ratio at 200 K is 1.5 
times that at 120 K. 

Helium Line Length 
ha = U5 m 

0.01 m < Ljj e < 2.0 m 

If the helium storage tank were 
mounted directly on the helium 
check valves, the fuels ide pre- 
burner would be unstable. Longer 
helium line lengths lead to more 
stable operation. Above 1.5 m, 
increasing the helium line length 
had only a small effect on oscil- 
lation frequency or amplitude 
ratio. 
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Table 4.3: (continued) 


ENGINE PARAMETER 
BASE VALUE 

RANGE OF VALUES USED IN 
TRNCHG 

EFFECT ON PREBURNER OPERATION 

Helium Line Diameter 

D„ = 0.739 cm 
He 

0.173 cm < D„ < 1.029 cm 
He 

Large helium line diameters 
volumes of helium in the line) 
lead to higher chamber pressures 
and oscillations which have higher 
frequencies and damp more slowly. 
For the smallest helium line dia- 
meter, a full solution could not 
be generated because numerical in- 
stability was encountered. 









I 


» 


ft 


CHAPTER V 
CONCLUSIONS 

The objectives of this study were to evaluate the validity of 
the TSTR model for modelling a rocket engine combustion chamber, to 
predict the amplitude and frequency of the SSME fuelside preburaer 
shutdown chug and to determine the processes and engine parameters 
that give rise to the chug. 

To meet the first two objectives the fuelside preburner was 
modelled as a heterogeneous TSTR combustion chamber, a variable mass 
flow rate oxidizer feed system, a constant mass flow rate fuel feed 
system and an exit turbine. A computer program, TRNCHG, was written 
to integrate the resulting differential equations numerically. For 
full power level simulations TRNCHG predicted stability and executed 
without difficulty. For SSME engine shutdown conditions, chemical 
kinetics were assumed infinitely fast and TRNCHG predicted quickly 
damped chamber pressure and oxidizer mass flow rate oscillations. 

To meet the third objective, engine parameters which had been 
shown in literature to influence the stability of liquid rocket 
engine combustion were varied. This procedure was quite successful, 
showing that small changes in some engine parameters resulted in 
large changes in prebumer operation and pointing out two engine 
parameters that have a larger role in the SSME fuelside preburner 
chug than predicted by Lira's (1986) linearized model. 
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Three questions are answered in this chapter: which engine 
parameters are the key parameters in the SSME fuelside preburner 
chug? how could the preburner model be extended or modified? and how 
well does the fuelside preburner model work? 

5.1 KEY ENGINE PARAMETERS AND THEIR EFFECT ON PREBURNER STABILITY 

As stated in Chapter IV, it is probably a combination of 
factors that result in the fuelside preburner chug. The most 
important factors and their influence on the prebumer chug, as 
predicted by the program TRNCHG, are given below. 

1. Mean droplet diameter. Mean droplet diameter is a function 
of injector geometry, chamber pressure and the ratio of the 
oxidizer and fuel mass flow rates. As mean droplet diameter is 
increased, pressure oscillations in the combustion chamber damp 
more slowly. This is due to an increase in the combustion delay 
time (or decrease in the vaporization rate) associated with 
large droplets. For a sufficiently large droplet ‘diameter, 
reaction blows off. 

2. Pressure downstream of the exit turbine (PHG). During engine 
shutdown, PHG falls from 2.41 MPa to 1.03 MPa (350 psia to 150 
psia). Pressure oscillations were shown to damp more slowly for 
lower values of PHG. 

3. Fuel temperature. As described in Chapter I, fuel cools the 


main combustion chamber and nozzle before being injected into 
the fuelside prebumer combustion chamber. During shutdown, the 
fuel mass flow rate to the combustion chamber is less than at 
full power level, but the main combustion chamber and nozzle 
are also cooler, so the temperature of the fuel injected into 
the fuelside prebumer may be lower during shutdown. Lower fuel 
temperatures result in a less stable fuelside prebumer. When 
the fuel temperature is reduced below 120 K, pressure 
oscillation frequency decreases rapidly and the oscillations 
become less damped. 

4. Fuel mass flow rate. Contrary to the predictions of Lim (1986), 
lower fuel mass flow rates were predicted to cause oscillations 
in chamber pressure to damp out more slowly. There are two 
reasons for this. First, for lower fuel mass flow rates there 
is less cold fuel to heat to the chamber temperature and more 
energy is available (from combustion) to sustain pressure 
oscillations. Second, the mean droplet diameter is a function 
of the fuel mass flow rate: as fuel mass flow rate decreases, 
the mean droplet diameter increases. As shown in figure 4.10 
(page 75), as mean droplet diameter is increased, preburner 
operation becomes less stable. 

5. Helium temperature. Changing the Helium temperature did not 
substantially change the amplitude of pressure oscillations, 
but did influence their frequency. Frequencies of pressure 


oscillations for helium temperatures of 120 K and 300 K were 
111 Hz and 175 Hz, respectively. 

6. Helium purge line geometry. The length and diameter of the pipe 
connecting the helium storage tank to the helium check valve 
both have a strong influence on the prebumer's operation 
during shutdown. As predicted by Summerfield (1951), as the 
length of the helium line was increased, the prebumer became 
more stable. Above a helium purge line length of about 1.0m, a 
large change in pipe length is necessary to accomplish a small 
change in frequency or the rate at which oscillations damp out. 
It was not possible to obtain solutions from TRNCHG for small 
helium purge line diameters. In general, small helium purge 
line diameters lead to reduced chamber pressures and greater 
stability. For a given helium line length, the volume of helium 
between the helium storage tank and the combustion chamber is 
important in determining whether fuelside prebumer will be 
stable 

Comparing the results generated by the non-linear model used in 
this study to the results generated by the linear model used by Lim 
(1986) to model the SSME fuelside preburner, the non-linear model 
showed the prebumer's stability to be sensitive to engine 
parameters to which the linearized model predicted the prebumer 
would be insensitive. One instance of the greater sensitivity of the 
non-linear model is that the inclusion of chemical kinetics in the 
calculation of reaction rates at full power level resulted in more 


slowly damped pressure oscillations than for chemical kinetics being 
considered infinitely fast (figure 4.28, page 100). Lim assumed that 
chemical kinetics were infinitely fast. Reaction rate is a 
non-linear function of chamber temperature and the mole number of 
the species in the combustion chamber. A second instance in which 
the predictions of the non-linear and linearized models differ is in 
the effect that fuel mass flow rate has on chamber pressure 
oscillations. Lim stated that chug could not be eliminated by 
changing the mass flow rate. As shown in figure 4.18 (page 86), the 
non-linear model predicts a pronounced effect of fuel mass flow rate 
on the rate at which oscillations are damped. Chamber temperature, 
droplet vaporization rate and mean droplet diameter are all 
functions of the fuel mass flow rate. 

Using a linearized model, Lim was able to predict the effect 
that a number of engine parameters had on the stability of the SSME 
fuels ide prebumer. The results of the linearized model were 
concise, required a relatively small amount of computer time and in 
most cases proved accurate. The non-linear model used in the current 
study was able to show that chemical kinetics and fuel mass flow 
rate play a greater role in determining whether or not the preburner 
undergoes stable operation than predicted by the linearized model. 
These predictions show that a non-linear model is necessary in a 
detailed study of chug and offset the greater computer time 
necessary to predict how the prebumer operates. 
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5.2 POSSIBLE EXTENSIONS OF THIS RESEARCH 


There are number of ways the fuels ide prebumer model could be 
changed to increase its accuracy and to enable it to be used to 
model different rocket engine phenomena. 

To improve the accuracy of the fuelside preburner model, the 
oxidizer feed system could be more faithfully represented, adding 
pipes and an oxidizer manifold. A finer grid of pipes would also 
increase accuracy, although for most cases this is not necessary. 
Within the oxidizer feed system the compressibility of LOX could be 
taken into account. 

In the combustion chamber, a droplet distribution, which, like 
Webber's (1972) could contain different sized droplets whose size is 
followed over time, could be added. Given temperature data for the 
combustion chamber walls, the rate of heat transfer from the 
combustion chamber walls to the chamber contents could be more 
accurately calculated. 

Finally, if it were desired to study spatial variations in 
chamber pressure, temperature and composition, it would be possible 
to use multiple TSTR's to model the combustion chamber. TSTR's could 
be distributed longitudinally and radially throughout the combustion 
chamber. The exit mass flow rate from the face of one TSTR would be 
an inlet mass flow rate through the face of an adjacent TSTR. Such a 
scheme could prevent the mixture ratio from being overest ima ted 


severely in sections of the combustion chamber that contain many 
liquid droplets. 

It is reiterated here that the above alterations of the TSTR 
model were not necessary in the study of the SSME fuelside preburner 
chug and are presented as ways in which the fuelside preburner model 
can be extended to be used in different applications. 

5.3 EVALUATION OF THE FUELSIDE PREBURNER MODEL 

The TSTR model of the fuelside prebumer combustion chamber is 
simple, yet has been shown capable of predicting low frequency 
combustion instability. With the proper assumptions, the model can 
be used over a wide domain of operating conditions and engine 
parameters . 

Among the strengths of the fuelside preburner model are, 

1. the frequency of pressure oscillations in the preburner is 
predicted to within 6% of typical measured frequencies, 

2. the effect that various engine parameters and combinations 
of these parameters has on preburner operation was shown, 

3. the presence of liquid in the combustion chamber is included, 

A. chemical kinetics may be included for high temperatures and 

pressures and 

5. the fuel prebumer is predicted to maintain stable operation at 
full power level. 


The preburner model's weaknesses are 

1. chemical kinetics are not modelled during the SSME shutdown, 

2. solutions cannot be generated when the diameter of the pipe 
connecting the helium storage tank to the helium check valve is 
too small, 

3. accurate initial conditions are required and 

4. long CPU times are required to generate solutions when large 
derivatives are encountered (e.g., when backflow of combustion 
chamber contents into the oxidizer feed system occurs). 

Despite the model's limitations, the objectives of this 

research were basically met. The TSTR model was able to predict 
transient conditions in a rocket engine prebumer. The model 
predicted that, during the SSME shutdown the fuelside preburner 
would not undergo unstable operation, yet when key engine parameters 
were changed by a small amount instability resulted. It is also 
possible, as a worst case, that the inclusion of chemical kinetics 
in the model for shutdown would result in unstable preburner 
operation. The frequency of the SSME preburner chug was accurately 
predicted, but the amplitude was not, since sustained oscillations 
were not predicted by TRNCHG. Finally, the model not only showed the 
influence that engine parameters and processes have on chug, but 
also pointed out, contrary to the predictions of a linearized model, 
that fuel mass flow rate and chemical kinetics are important. 

The program, TRNCHG, was written for the SSME fuelside 
preburner. However, by altering the oxidizer feed system 
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calculations and changing the way that exit mass flow rate is 
calculated, TRNCHG could be applied to many other rocket engine 
combustion chambers to study low frequency combustion instability. 
Using such a program has been shown valuable in an in-depth study of 
chug. 
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APPENDIXES 


APPENDIX A 


COMBUSTION CHAMBER AND OXIDIZER FEED 


SYSTEM DERIVATIVES 


A. 1 DERIVATION OF THE TSTR GOVERNING EQUATIONS 


The conservation of mass for the gas and the liquid phases is 
fully derived in Chapter III. The conservation of species and energy 
for a TSTR and the ideal gas law as used in the present study are 
derived below. 

The Conservation of Species 


Equation [3.6] can be written 


d(o i &) 


NSTRM 


dt 


E fm. o. . ^ - m a. + T. 
j =1 \ JV ij/ gv 1 1 


[A. 1 ] 


for i=l,2,3 ...NS. The subscript v denotes a volumetric quantity. 
The derivative in [A.l] is expanded, giving 


NSTRM, ^ ^ 


„ d °i d Pc / •* * \ . 

P — - + o. — = E { m. a. .) - m a. + T. 

^ dt 1 dt j=1 V jv ij; gv i iv 


From equation [3.4], 



NSTRM 



[A. 2] 


[A. 3] 


[A. 3] is substituted into [A. 2] and the result is simplified, 
yielding 
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da. 

l _ 

1 

NSTRM 

E 

dt 

Pc 

II 


IV 


tA. 4] 


For any variable, x, 


d(ln x) 1 dx 

dt x dt 


[A. 5] 


[A. 5] is applied to [A. 4], giving the conservation of species for 
the TSTR in its final form: 


d(ln o^) 1 


dt 


’iPc 


NSTRM 

E [m* (a*. - o.)] + 

j=1 Ijv ij i J 


IV 


i=l,2,3 ... NS 


[A. 6] 


Conservation of Energy 


Conservation of energy is given by equation [3.7] as 


dE 

dt 


NS NSTRM 


NS 


( .* * * \ /. \ . 
ra.o..h..)- Elm o.h. )- Q 

j u uJ i=1 \ s i iJ 


[A. 7] 


assuming that changes in the kinetic and potential energies of the 
chamber contents are negligible. 


dE 

dt 



[A. 8] 


where E is the total energy (Joules) of the chamber contents, U is 
the internal energy (Joules) of the chamber contents and u^ is the 


specific molar internal energy (J/kgmole) of species i. Expanding 
the derivative, 


[A. 9] 


dE 

NS NS 

do. 

NS 

du. 

— = 

* — V 2 u.o. + P V 2 

1. 

u. 

+ p V 2 o. 

i 

dt 

-ill . . 

dt i=l i=l 

dt 1 

i-1 1 

dt 


Assuming that the contents of the chamber behave as a mixture of 
perfect gases and that their specific heats at constant pressure, 
c ^ ( J/kgmole *K) do not vary. 


du. d „ dT ■ x 
— - = — (h. • R T ] = — - (c . - R ) 
dt dt ' 1 a* \ P 1 SJ 


dT 

c 

dt 


[A. 10] 


Using [A. 10], [A. 9] becomes 
dE d p c NS 


NS do. 


dt dt 


V 2 o.u. +OV 2 — - u. + 
i 1 r'c l 

1=1 1=1 dt 


dT NS 

P V — - 2 o.(c . - R ) 

Hc dt i=l 1 P 1 * 


[A. 11] 


[A. 11] is substituted into [A. 7] and [A. 5] is applied, giving 
•NS NSTRM NS 5 NS 


d( In T 


dt 


v i-imo iNoiKn no 

. ' . * A * . d P c 

— = 2 2 m. o..h..-2m o.h. - 2 o.u. 

L i=i j-i 1J 1J i=i 1 1 dt i=i 1 1 


NS do NS 

Cl 2 — u. - Q/V / T P 2 o. ( 

^ i=l dt 1 J [ c i=l 1 


c . - R ) 
Pi g 


[A. 12] 
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The Ideal Gas Law 


For the combustion chamber, the ideal gas law can be written 


P = D a R T 
c r'c m g i 


g c 


NS 


where o = E a. . Alternatively, 


m i=l 1 


dP d 

c ^ „ c 


- + D ° R 

m g c r'c m g 


= a R T 

dt m g c dt 


dT do 

— - + p R T — 
dt Hc 8 c dt 


[A. 13] 


Equation [A. 5] is applied to [A. 13], giving 


d(ln P ) oRT dfi floR dT DR T do 
c __ m g c 'c t re m g c | ^c g c m [A. 14] 


dt 


dt 


P dt 
c 


P dt 
c 


or. 


d(ln P ) 1 dn . 1 dT 1 do 

c re + c , m 


dt 


Pc dt 


T dt 
c 


o dt 
m 


[A. 15] 


A. 2 DERIVATION OF THE REACTION RATE FROM ELEMENTARY KINETICS 


The volumetric forward and reverse reaction rates for an 
elementary reaction, j, can be written (Pratt and Wormeck, 1976) 



RV 


jv 


M ! , NS , 

Aj r; 1 expf-r /T c ) (fto/m , (p c <,.)1i [A. 16) 

1=1 

n NS II 

T^iexpl-T" /T c ) , (fto )*i| [A. 17) 

1=1 
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where A., N. and T . are empirically determined coefficients, T is 
J J aj c 

the combustion chamber temperature, is the stoichiometric 

coefficient of species i in reaction j , u is the stoichiometric 

m 

coefficient for a third body, and the superscripts ' and " denote 
forward and reverse reaction, respectively. The overall rate of 
production of species i by a chain of NRXN elementary reactions is — 


NRXN 

R. = E (v" .R" - v\ .R! ) 

j=l 


IV 


1J JV 


1J JV 


[A. 18] 


The following elementary reactions made up the overall 
conversion of and 0 2 to water vapor in the TSTR model of the SSME 
fuels ide preburner: 


OH + H 2 
0 2 + H 

0 + h 2 
0 + h 2 o 

H + H + M 
0 + 0 + M 
0 + H + M 
H + OH + M 

Because the atoms 0 and H and the 


< — -> OH + OH 

<---> H 2 0 + H 
<---> OH + 0 
<---> OH + H 
<---> OH + OH 
<---> H 2 + M 
<---> 0 2 + M 

< — -> OH + M 

< — -> H 2 0 + M 

radical OH are so reactive, they 


are consumed nearly as soon as they are produced, so their mole 


numbers were assumed to be constant at any given time, or 


^0 = = dq 0H 

dt dt dt 


= 0 


[A. 19] 


Using equation [A. 19], equation [A. 4] can be written for 0, H and 
OH, giving 


NSTRM 


f H ‘ ^ "jv (o Hj ' V + V, = 0 


NSTRM 


.* 


f _ = £ m. (a^ . - o A ) + R_ = 0 

0 jv Oj 0 Ov 


[A. 20] 


NSTRM * 

f OH ' V ( °0Hj ' °0H ) + E 0Hv * 0 


Since there are no radicals present in the inlet streams. 


f H " -'V s "jv 5 + "hv = 0 


.* 


f 0 * -V 1 "jv 5 + R 0v ’ 0 


• * 


f 0H = "°0H (S "jv 5 + R 0Hv * 0 


[A. 21] 


At each time, values for o„ , a n , a„ n , T , P and are known 

H 0 H 0 c c c 

(these are the combustion chamber variables which are being 

integrated), so equations [A. 21] constitute three equations in the 

three unknown mole numbers o„, a n and o™. 

H U Un 


The roots of equations [A. 21] were found using a Newton 
iteration. The functions, f^, were expanded in a Taylor series 
about the assumed solution and second order and higher terms were 
truncated (George, 1982), giving 


f. 

1 


= f i,old 


NRAD 

+ E 
k=l 


3 f 


i.old 



i=l ,NRAD [A. 22] 


where NRAD is the number of species whose mole number is assumed 
constant at each time step (Nrad=3 for the SSME TSTR model). If the 
correct values of the radical mole numbers are chosen, f.= 0 for 

l 

i=l,NRAD. [A. 22] can then be written 
NRAD 0 f 

E Ao, - -f . . , i=l ,NRAD 

, i P k i,old 
k*l Go k 

or. 


dydon <3ydc-o df E /do QE 



/ \ 
" f H,old 

df 0 /do E df 0 /do Q df 0 /do 0n 

o 

<1 

> = i 

“ f O,old 

^OH ^°H ^ f 0H /( ^ a 0 ^ f 0H /C ^ O 0H_ 

a °oh 


' f OH, old 


Equation [A. 23] is solved for the correction variables 
i=l,NRAD. New estimates for o^, i=l,NRAD, are 


a. -a. + Ao. 

i,new i,old i 


[A. 24] 


The results of [A. 24] are successively substituted into [A. 23] until 


the condition 


o. < e i=l ,NRAD [A. 25] 

1 ss 

is met. e ss was set at 1.0x10 ^ in the program TRNCHG. 

When condition [A. 25] was met, reaction rate for H^O was 
calculated using the known values for mole numbers, temperature, 
pressure and density, the steady state values for the radical mole 
numbers and equations [A. 16], [A. 17] and [A. 18]. By conservation of 
atoms , 

’HijV = " ^HjOv 

V * - °- 5 Vov ' A - 26 l 

In all observed cases, the approximate values for Rg v and R Q ^ were 
close (within an estimated 5%) of the values calculated directly 
using [A. 16], [A. 17] and [A. 18]. 

A. 3 DERIVATION OF THE OXIDIZER FEED SYSTEM EQUATIONS 


In integral form, the conservation of linear momentum for a 
control volume is (White, 1979) 


EF = T* ( JJJ P vdV ) + JJpv(vn) dA [A. 27] 


As shown in figure 3.6, using equations [3.15] and [3.16], 


EF = (P - P 2 )A - fLPv 2 A 
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For a pipe containing a liquid/ vapor interface, [A. 27] can be 
written 

d 

(Pi - P 2 M - fLpv 2 A = — (p^v + Pg V g V ) + P^A - P g v 2 A [A. 29] 
Note that 



dt dt 


[A. 30] 


Divide by A and group terms for 


dv d A d D 

P - P - fLpv 2 = (P^l + P ' 1 ) ~ + V + 1 v 

8 g dt dt 8 dt 


[A. 31] 


For each pipe node, k, equation [3.20] is written 


NPN NPN 

P k Z A n (dv n /dt) + (d^/dt) E A n v n = 0 [A. 32] 

n=l n=l 


Equations [A. 31] and [A. 32] constitute NN+NP linear equations in 
the NN+NP unknowns: dv n /dt, n=l,2,3 ...NP, and cp^/dt, 
k=l,2,3, ...NN. [A. 31] and [A. 32] can be expressed as the matrix 
equation [A] {X} = (B) where matrix [A] contains the coefficients of 
the unknowns, (X) contains the values of the unknowns and (B) 
contains " ^2 ” ^P v2 for rows 1,2,3 ...NP and zero for rows 
NP+l,NP+2, . . . NP+NN. In the program written for this study, (X) was 
found via a call to IMSL subroutine LEQT1F. 


APPENDIX B 


FLOWCHARTS, PROGRAM LISTING AND SAMPLE INPUT 
AND OUTPUT FILES FOR PROGRAM TRNCHG 
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TRNCHG : MAIN (AUGUST, 1987) 


C 
C 

C TRNCHG WAS WRITTEN TO PREDICT THE TRANSIENT BEHAVIOR OF THE 
C SPACE SHUTTLE MAIN ENGINE FUELS IDE PREBURNER DURING THE SHUT- 
C DOWN HELIUM PURGE. IT IS ALSO CAPABLE OF PREDICTING ENGINE BE- 
C HAVIOR DURING OPERATION AT FULL POWER LEVEL. 

C 

C INITIAL PROGRAMMING FOR TRNCHG WAS BY DR. P. GEORGE II, CURRENT - 
C LY AT BATTELLE LAB, COLUMBUS, OH. SUBSEQUENT PROGRAMMING WAS BY 
C T. BARTRAND, MASTER'S CANDIDATE, U. OF TENNESSEE, KNOXVILLE. 

C 

C THE PREBURNER COMBUSTION CHAMBER WAS MODELLED AS A STIRRED TANK 
C REACTOR. FUEL MASS FLOW RATE WAS TAKEN AS CONSTANT. OXIDIZER EN- 
C TERS THE COMBUSTION CHAMBER VIA A MULTI-PIPE MULTI-NODE FEED 
C SYSTEM WITH A CONSTANT PRESSURE AT THE FURTHEST UPSTREAM NODE. 

C 

C THE TRNCHG SYSTEM VARIABLES, WHICH ARE INTEGRATED OVER TIME, 

C FALL INTO TWO CATEGORIES: THE' CHAMBER VARIABLES AND THE OXI- 
C DIZER FEED SYSTEM VARIABLES. 

C 

C CHAMBER VARIABLES INCLUDE: 

C S2(I) : MOLE NUNMBER OF SPECIES I (KGMOLES OF I /KG) 

C TK : TEMPERATURE (K) 

C PA : PRESSURE (PA) 

C DLIQ : LIQUID DENSITY IN THE COMBUSTION CHAMBER (KG/M**3) 

C 

C FEED SYSTEM VARIABLES INCLUDE: 

C V(J) : VELOCITY IN PIPES J (M/S) 

C FVI(J) : POSITION OF THE LI QUID/ VAPOR INTERFACE IN PIPE J 

C (% OF THE PIPE LENGTH) 

C RHO(K) : DENSITY AT NODE K (KG/M**3) 

C FVI2 : POSITION OF THE BACKFLOW INTERFACE IN THE LAST PIPE. 

C (% OF THE PIPE LENGTH) 

C 

C EXECUTION OF TRNCHG INVOLVES 11 SUBROUTINES (WRITTEN FOR THE 
C PROGRAM) AND THE USE OF TWO IMSL SUBROUTINES: DGEAR, FOR NUMER- 
C ICAL INTEGRATION OF DIFFERENTIAL EQUATIONS AND LEQT1F, FOR SOLU- 
C TION OF MATRIX EQUATION A X = B FOR THE VECTOR X. A LISTING OF 
C THE SUBROUTINES USED IN THE EXECUTION OF TRNCHG AND A BRIEF 
C EXPLANATION OF THE FUNCTION THAT EACH PERFORMS FOLLOWS. 

C 

C COMMONSS : COMMON BLOCK INCLUDED IN MOST OF THE SUBROUTINES 
C CALLED IN THE EXECUTION OF TRNCHG. 

C YCAL : USED TO CALL THE INTEGRATION ROUTINE AND CALL SUBROUTINE 
C OUTPT . 

C NPT: READS MOST OF THE NECESSARY DATA FOR EXECUTION. INITIAL - 
C IZES VARIABLES, CALCULATES THERMODYNAMIC PROPERTIES AND 

C FORMS REACTION RATE EXPRESSIONS. 

C DTSTR : RUNGE-KUTTA INTEGRATION SUBROUTINE TO BE USED IN THE 
C ABSENCE OF IMSL SUBROUTINE DGEAR. TO SELECT DTSTR THE VARI- 
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c 

c 

c 
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ABLE NDGEAR (NAMELIST ICONS) IS SET EQUAL TO 1 . TO USE DGEAR 
LET NDGEAR = 0. 

CSLOPE : CALCULATES THE DERIVATIVES WITH RESPECT TO TIME FOR 
THE SYSTEM VARIABLES ON A CALL FROM THE INTEGRATING ROUTINE. 
CALLS SUBROUTINES COMBRT , PIPES AND THERPS. 

THERPS: CALCULATES ENTHALPY, INTERNAL ENERGY AND SPECIFIC HEAT 
AT CONSTANT PRESSURE FOR GASEOUS CONSTITUENTS WITHIN THE 
COMBUSTION CHAMBER. 

COMBRT: CALCULATES THE RATE OF PRODUCTION OF THE MAJOR SPECIES 
IN THE COMBUSTION CHAMBER DUE TO VAPORIZATION AND COMBUSTION. 
COMBRT CALLS BOTH RADSS AND THERPS. 

PIPES: CALCULATES THE DERIVATIVES OF THE PIPING SYSTEM VARI- 
ABLES ON A CALL FROM CSLOPE. ACCESSES IMSL SUBROUTINE LEQT1F. 
INPUTS TO PIPES ARE FOUND IN NAMELIST INPIP AND COMMON BLOCK 
CPIPE. 

RADSS: CALCULATES THE RATE OF PRODUCTION OF THE MAJOR SPECIES IN 
THE COMBUSTION CHAMBER DUE TO COMBUSTION. ASSUMES THE MOLE 
NUMBERS OF THE VERY REACTIVE SPECIES (ATOMS AND RADICALS) ARE 
CONSTANT AT ANY TIME. RADSS ACCESSES IMSL SUBROUTINE LEQT1F 
AND SUBROUTINE RATES. 

RATES: CALCULATES THE RATE OF PRODUCTION OF RADICALS AND ATOMS 
AND THE DERIVATIVE OF THESE RATES. RATE FOR A GIVEN SPECIES 
IS FOUND BY SUMMING THE CONTRIBUTIONS TO RATE BY EACH ELEMEN- 
TARY REACTION. 

OUTPT: OUTPUTS VALUES TO THE TEXT FILE TSTR.OUT AND THE PLOT 
FILES TSTR.PLT1 AND TSTR.PLT2 ON RETURN FROM THE INTEGRATION 
SUBROUTINE OR IF THE EXECUTION MUST BE STOPPED 

OUTPUT IS CONTROLLED BY THE PARAMETER NDEBUG. FOR NDEBUG(I) EQUAL 

TO 1 THE OUTPUT OF OUTPUT GROUP I IS PRINTED IN THE APPROPRIATE 

FILE. OTHERWISE NDEBUG(I) SHOULD EQUAL 0. 


OUTPUT GROUPS ARE: 


NDEBUG ( 1 ) ---> 

NDEBUG(2) ---> 
NDEBUG(3) ---> 
NDEBUG(4) ---> 

NDEBUG(5) ---> 
NDEBUG(6) ---> 

NDEBUG ( 7 ) ---> 

NDEBUG(8) ---> 

NDEBUG(9) ---> 

NDEBUG ( 10) --> 


PROBLEM PARAMETERS FROM NAMELISTS PARAM, ICONS, 
SSRXN AND PIPES. 

THERMODYNAMIC PROPERTY DATA. 

REACTION RATE DATA. 

INITIAL CONDITIONS OF THE INLET STREAMS 
AND IN THE COMBUSTION CHAMBER. 

INPUTS TO THE INTEGRATING ROUTINE. 

OUTPUTS FROM THE INTEGRATING ROUTINE TO 
THE OUTPUT TEXT FILE (TSTR.OUT) 

OUTPUTS FROM THE INTEGRATING ROUTINE TO 
THE OUTPUT PLOT FILE (TSTR.PLT) 

TERMS USED IN THE CALCULATION OF DERIVATIVES 
TIVES IN SUBROUTINE CSLOPE. 

DERIVATIVE ESTIMATES AT THE END OF 
SUBROUTINE CSLOPE. 

OUTPUTS VALUES OF TEMP, PRESS, AND MOLE NO. TO 


» 
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NDEBUG(ll) 

NDEBUG(12) 

NDEBUG( 13) 
NDEBUG( 14) 

NDEBUG(15) 


THE SCREEN IN CSLOPE (FOR INTERACTIVE USE) 

--> CONVERGENCE SUMMARY FOR SUBROUTINE RADSS 

--> SUMMARY OF THE DERIVATIVES CALCULATED IN 
SUBROUTINE PIPES. 

--> OUTPUTS FROM SUBROUTINE RATES, INCLUDING 
A LISTING OF THE RXNS WHOSE REACTION RATES 
ARE SMALL AND SET EQUAL TO ZERO. 

--> OUTPUTS VALUES OF TIME, STEP, TEMPERATURE 
AND MOLE NUMBERS OF THE MAJOR SPECIES TO 
THE SCREEN ON RETURN TO YCAL FROM THE INTEG- 
RATION SUBROUTINE. 

--> SUMMARY OF DERIVATIVE ESTIMATES, OLD AND 
NEW SYSTEM VARIABLE VALUES FROM THE RUNGE - 
KUTTA INTEGRATION SUBROUTINE DTSTR. 


THREE SETS OF COMBUSTION CHAMBER CONSTITUENTS ARE 
CONSIDERED IN TRNCHG: 

(1) INLET STREAM CONSTITUENTS (REACTANTS) WHICH ARE 
PRESENT IN SIGNIFICANT AMOUNTS IN THE COMBUSTION 
CHAMBER, 

(2) REACTION PRODUCTS WHICH ARE PRESENT IN SIGNIFICANT 
AMOUNTS IN THE COMBUSTION CHAMBER AND 

(3) RADICALS AND ATOMS, WHICH ARE CONSIDERED TO HAVE 
SMALL CONCENTRATIONS RELATIVE TO REACTANTS AND 
PRODUCTS AND WHOSE MOLE NUMBER IS TAKEN CONSTANT 
AT ANY INSTANT OF TIME. 


TRNCHG IS DOUBLE PRECISION FOR ALL REAL VARIBLES . THE MAXI- 
MUM NUMBER OF SYSTEM VARIBLES IS 50, BUT CAN BE EASILY 
INCREASED BY INCREASING THE DIMENSIONS OF ARRAYS Y, YPRIME, 
WK AND IWK. 

***** TRNCHG USES SI FOR ALL CALCULATIONS ***** 

INCLUDE 'INIT. FOR/LIST ' 

DIMENSION S2P( 10) ,Y(50) ,IWK(50) ,WK( 9000) 

INCLUDE 'COMMONSS. FOR/ LIST’ 


NAMELIST 

* /RUNNO/ 

* /PARAM/ 

* 

* 

* /INITC/ 

* /INDX1/ 

* /ICONS/ 


NMON , NDAY , NYR , NRUN , 

VOL , Q , S2MIN , EMS , TKS , TMSD , RGMX , RHOGAS , RHOLIQ , 
RHOH2 , AH2IN , A02IN , FVI2MN , DLIQMN , DMMIN , DMC , 
CFTP , PHG , XPER , TMPER , CVR , 

XMASS,TK, PA,DM, 

NDEBUG , MCON , KINET , 

METH , MITER , NDGEAR , IER , STEP , TMI , TMPRNT , TMEND , 


140 


onoon on ooo 


* INDEX, EPSI, 

* /SSKRT/ EPS ISS , RLX , RLX2 , ITMAXS , IDGT , JER 

* /INPIP/ V,FVI,RHO, FVI2.PL, AREA, RFLIQ,RFGAS,RFLHO, 

* RFGHO ,RFLV,RFGV,PS, CDIN , SM1MIN , 

* /INDX2/ N1 ,N2,NP,NN,NHE0,NFVI 

• C 

DATA RGAS/ 8314.4 /,JJ/ 9 / ,TENLN/ 2.302585093 / 

C 

OPEN (UNIT=20, FILE*' [BART. TSTR. WORK]TSTR. IN' .READ ONLY, 

* STATUS*' OLD 1 ) 

OPEN ( UNIT=3 0 , FILE* ' [ BART. TSTR. WORK ]TSTR. OUT’ , STATUS* ' OLD ' ) 
OPEN (UNIT=31 .FILE* ' [ BART. TSTR. WORK] TSTR. PLT1' .STATUS* 'NEW' ) 
OPEN ( UNIT=32, FILE* ' [ BART . TSTR . WORK ] TSTR . PLT2 ' , STATUS* ' NEW ' ) 

--- READ AND WRITE PROBLEM PARAMETERS FROM NAMELISTS. 

READ (20.RUNNO) 

READ (20.PARAM) 

READ (20.INITC) 

READ (20.INDX1) 

READ (20, ICONS) 

READ (20.SSRRT) 

READ (20, INPIP) 

READ (20.INDX2) 

WRITE (30,1000) NMON , NDAY , NYR , NRUN 
IF (NDEBUG(l)) THEN 

WRITE (30,1100) VOL , Q , S2MIN , XMASS , DM , DMMIN , CFTP , PHG , 

* NDEBUG 

WRITE (30,1150) DMC , AH2 IN , A02 IN , MCON , KINET 

WRITE (30,1200) METH, MITER, NDGEAR.IER, EPSI, STEP, TMI, 

* TMPRNT , TMEND , TMSD 

WRITE (30,1300) EPS ISS, RLX, ITMAXS, IDGT, JER 

WRITE (30,1400) NP,NN,NHEO,NFVI,RFLIQ,RFGAS,RFLHO, 

* RFGHO , RFLV , RFGV , CDIN , PS , FVI 2 

WRITE (30,1500) (I,N1(I) ,N2(I) ,AREA(I) ,PL( I) ,1*1 ,NP) 

WRITE (30,1600) (I,V(I) ,FVI( I) ,1=1 ,NP) 

WRITE (30,1700) (I,RHO(I),I=l,NN) 

END IF 

S2MLG = DLOG(S2MIN) 

PAMAX * 1.0D+03 * PA 

— ARATIO IS THE RATIO OF THE TOTAL AREA OF THE OXYGEN 
INJECTORS TO THE AREA OF THE LAST PIPE IN THE OXIDIZER 
FEED SYSTEM. 

ARATIO = 264.0 * A02IN/AREA(NP) 

ARATSQ = ARATIO*ARATIO 
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--- CALL NPT TO ESTABLISH INITIAL CONDITIONS FOR THE REACTANT 
STREAMS AND COMBUSTION CHAMBER. NPT INDEXES SPECIES AND 
FILLS S2 AND SI ARRAYS AND INITIALIZES TEMPERATURE AND 
PRESSURE. 

CALL NPT 

--- DEFINE THE IS IDE ARRAY. ISIDE(KK.LL) TELLS THE SIDE OF 
RXN LL ON WHICHy SPECIE KK IS FOUND. IF KK ON THE LHS, 
ISIDE IS NEGATIVE. IF KK IS ON THE RHS, KK IS POSITIVE. 

IF KK APPEARS TWICE, ISIDE-2. IF KK APPEARS ONCE, 

ISIDE=1. FOR KK NOT PRESENT, ISIDE=0. 

DO 120 LL=1 , JJ 
DO 110 KK=1,NSRATE 
DO 100 MM=1 , 4 

IF(KK.EQ.ID(MM,LL)) THEN 

IF(MM.LE.2.AND.ID(1,LL).EQ.ID(2,LL)) THEN 
ISIDE(KK.LL) = -2 
ELSE IF(MM.LE. 2) THEN 
ISIDE(KK,LL) = -1 

ELSE IF(ID(3,LL).EQ.ID(4,LL)) THEN 
ISIDE(KK,LL) = 2 
ELSE 

ISIDE(KK.LL) = 1 
END IF 
GO TO 110 
END IF 
100 CONTINUE 

ISIDE(KK.LL) = 0 
110 CONTINUE 
120 CONTINUE 


--- WRITE THERMODYNAMIC AND REACTION RATE DATA TO THE OUTPUT 
TEXT FILE. 


IF(NDEBUG(2)) THEN 

WRITE (30,1900) (I,ANAM(I),SMW(I),I=1,NS) 

WRITE (30,2000) (I,ANAM(l),SMW(I),I=NS+l,NSRATE) 
WRITE (30,2100) (ANAM(I),(Z( J),J=(I-1)*14+1,I*14), 
I-l.NSRATE) 


END IF 

IF(NDEBUG(3)) THEN 

WRITE (30,2200) (J,(ID(K,J),K=1,4), J=1,JJ) 

WRITE (30,2300) ( J, (ISIDE(K, J) ,K=1 ,NSRATE) , J=1 , JJ) 
WRITE (30,2400) 

DO 125 J=1,JJ 

TM1 = BX( J)/TENLN 
TM2 = BX2(J)/TENLN 

WRITE (30,2450) J,MODR( J) ,TM1 ,TEN( J),TACT( J),TM2, 


TEN2 ( J ) , TACT2 ( J ) 


* 


C 

C 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 


c 


c 

c 


125 CONTINUE 
END IF 

---PUT MOLE NUMBER WITHIN THE PROBLEM DOMAIN 

(S2MIN< S2 < 1/SMW), SUMMARIZE INPUTS AND WRITE INTIIAL 
CONDITIONS TO THE PLOT AND TEXT OUTPUT FILES. 


SM = 0.0 
DO 130 1=1, NS 

IF(S2(I) .LT.S2MIN) S2(I) = S2MIN 
SM = SM + S2(I) 

130 CONTINUE 

--- DEFINE USEFUL CONSTANTS 

TKINV = 1.0/TK 

SMINV = 1.0/SM 

RGASIN = 1.0/RGAS 

OMV = EMS( l)/VOL 

FMV = EMS(2)/VOL 

PMV - EMS(3)/VOL 

EMV = OMV + FMV + PMV 

RHOCC = PA*SMINV*RGASIN*TKINV 

WRITE INITIAL PIPE AND COMBUSTION CHAMBER CONDITIONS TO 

THE TEXT FILE AND THE PLOT FILE. 

IF(NDEBUGU)) THEN 

WRITE (30,2500) (I,ANAM(I),S2(I),ISTRM(I),S1(ISTRM(I), 

* I ) , AST( I ) , 1= 1 , NSRATE ) 

WRITE (30,2600) TK, PA, RHOCC, SM, VOL, (I,ANST(I) ,TKS(I) , 

* EMS(I) ,SF( I) ,HS(I) ,SH(I) ,1=1 ,3) 

END IF 

IF(NDEBUG(7)) THEN 
TMP = TMI*1 . 0D+03 
PAP = PA*1 . 0D-06 
IF(DM.LT.DMMIN) DM = DMMIN 
DMP = DM * 1.0D+06 
DLIQ = XMASS/VOL 
DO 140 1=1, NSRATE 
140 S2P(I) = S2(I)*1 . 0D+03 

WRITE (31,2700) 

WRITE (31,2800) TMP, (S2P(I) ,1=1 , NSRATE) 

WRITE (32,2750) 

WRITE (32,2850) TMP, TK, PAP, EMS(l), DLIQ, DMP, VR 
ENDIF 

SET PIPES INITIAL CONDITIONS AND INDICIES. NOTE THAT 
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THE LOG VARIABLE FORM IS NOT USED IN THE FORMULATION 
OF THE PIPING SYSTEM OF EQUATIONS. 

NN NO. OF NODES 

NP NO. OF PIPES 

NHEO NO. OF PIPES THAT MAKE UP THE HELIUM PURGE 
ORIFICE 

NPV INDEX FOR THE PIPE THAT ACTS AS THE HELIUM 
CHECK VALVE 

NFVI INDEX FOR THE PIPE IN WHICH THE LIQUID/VAPOUR 
INTERFACE IS LOCATED 

NS4 --> NS 5 INDICIES FOR PIPE FLUID VELOCITIES. 

NS6 --> NS7 INDICIES FOR PIPE FLUID/VAPOUR INTERFACES. 

NS 8 — > NS9 INDICIES FOR NODAL DENSITIES 
NS 10 INDEX OF THE SYSTEM VARIABLE THAT FOLLOWS THE 
PROGRESS OF A LIQUID/VAPOR INTERFACE THAT RE- 
SULTS FROM BACKFLOW FROM THE COMBUSTION CHAM- 
BER INTO THE OXIDIZER FEED SYSTEM. 

NV TOTAL NUMBER OF DEPENDENT (SYSTEM) VARIABLES 

WHICH MUST BE INTEGRATED. 

NODE 1 IS THE FURTHEST UPSTREAM NODE. PRESSURE FOR NODE 
1 IS TAKEN AS CONSTANT. NODE NN IS THE FURTHEST DOWN- 
STREAM NODE. IT IS THE JUNCTION OF THE PIPING SYSTEM WITH 
THE COMBUSTION CHAMBER. ONLY ONE UPSTREAM AND ONE DOWN- 
STREAM BOUNDARY ARE PROGRAMMED. 

PIPE 1 IS TAKEN AS THE FURTHEST UPSTREAM PIPE. PIPE NP IS 
THE FURTHEST DOWNSTREAM PIPE, CONNECTED TO THE COMBUSTION 
CHAMBER VIA NODE NN. 

NS4 = NS3 + 1 
NS5 = NS3 + NP 
NS 6 = NS5 + 1 
NS 7 = NS5 + NP 
NS 8 = NS7 + 1 
NS 9 = NS7 + NN 
NS10 = NS9 + 1 
NV = NS10 
NPV = NHEO + 1 

N AND NWK ARE THE DIMENSIONS OF THE SYSTEM VARIABLES USED 
IN THE INTEGRATION. 


N = NV 

IF(METH.EQ.l) NMETH = NV*13 
IF(METH.EQ.2) NMETH = NV*6 

IF ( MITER . EQ . 1 . OR . MITER . EQ . 2 ) NMITER = NV*(NV+1) 
IF(MITER. EQ. 3) NMITER = NV 

IF ( MITER . EQ . - 1 . OR . MITER . EQ . - 2 ) MITER=NV* ( 2*NLC+NUC+3 ) 
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IF ( MITER. EQ.O) NMITER = 1 
NWK = 4*NV + NMETH + NMITER 


--- YCAL CARRIES OUT THE INTEGRATION AND OUTPUTS RESULTS TO A 
PLOT FILE AND TEXT FILE ON RETURNS FROM THE INTEGRATION 
ROUTINE. 

CALL YCAL(Y,WK,IWK,N,NWK) 

C 

STOP 

C 

1000 FORMAT(//2X, ****** TSTR OUTPUT TEXT FILE *****'/, 

* 12X,I2,2X,I2,2X,I4/ , 12X, 'RUN NUMBER' ,2X, 13,//) 

1100 FORMAT(2X, ' ***** SOLUTION PARAMETERS *****'//, 

* 4X,'CC VOLUME = ',1PD9.2,' M**3'/, 

* 4X, 'CC HEAT TRANS. RATE = '.1PD12.5,' J/S'/, 

* 4X, 'MIN ALLOWABLE MOLE NO. = '.1PD12.5,' KGMOLE/KG' / , 

* 4X, 'INITIAL MASS OF LIQ IN THE COMB CHAMBER » \1PD12.5, 

* ' KG 7, 4X, 'INITIAL DROPLET DIAMETER = ',1PD12.5,' M 7, 

* 4X, 'MIN ALLOWABLE DROPLET DIAMETER = ’.1PD10.3,' M'/, 

* 4X , ' EXIT TURBINE FLOW RATE CONSTANT « '.1PD10.3/, 

* 4X, 'HOT GAS MANIFOLD PRESSURE = ',1PD10.3,' PA'/, 

* / 4X , ' NDEBUG (1 FOR PRINITING, 0 FOR SUPRESSING OUTPUT)' 

* / , (8X, 10(11, 2X))) 

1150 FORMAT (/4X, 'DROPLET DIAMETER COEFFICIENT = ',1PD11.4/, 

* 4X, 'AREA OF HYDROGEN INJECTOR ANULUS = ',D11.4/, 

* 4X, 'AREA OF THE OXIDIZER INJECTOR = ' ,1PD11.4/,4X, 

* 'MCON (=0 FOR CONSTANT MASS FLOW; =1 FOR PIPING SYSTEM) ' , 

* 11/ ,4X, 'KINET (=0 FOR NO KINETICS; =1 FOR CALL TO RADSS) ' 

* , 11 ) 

1200 FORMAT (/4X,' INTEGRATION PARAMETERS:'/, 

* 6X, 'MITER, METH, NDGEAR' ,2X,3(I1,3X)/,6X, 'IER, EPSI', 

* 2X,I3,2X,1PD9.2/,6X, 'INPUT STEP SIZE = ',1PD9.2,' SEC'/, 

* 6X, 'INITIAL TIME = ',D9.2,' SEC 7, 6X, 'OUTPUT INTERVAL', 

* ’ = ' ,D9.2, ' SEC 7, 6X, 'FINAL TIME = ',D9.2,' SEC'/, 

* 6X , ' TIME SHUTDOWN BEGINS = ’,D9.2) 

1300 FORMAT ( /4X, ' FOR SS RATE APPROXIMATION : 7 , 6X , 

* 'CONVERGENCE CRITERION (MAX % CHANGE) = ' , 1PD9.2/.6X, 

* 'RELAXATION PARAMETER = ' ,0PF10.5/ ,6X, 'MAX NO ITERATIONS', 

* ' = ' ,17/, 6X,' ERROR CHECK IDGT = ',12/, 

* 6X , ' IMSL ERROR FLAG JER = ’,13) 

1400 FORMAT (/4X, 'OXIDIZER FEED SYSTEM INPUTS:'/, 

* 6X, 'NO. OF PIPES', 2X, 12, 2X, 'NO. OF NODES' ,2X, 12/ , 

* 6X, ' NO. OF HELIUM ORIFICE PIPES’ ,2X, 12/, 

* 6X, ' PIPE IN WHICH THE INTERFACE IS INITIALLY LOCATED', 

* 2X,I2/,6X, 'RFLIQ FOR THE OX PIPES = ',1PD10.3/, 

* 6X, ' RFGAS FOR THE OX PIPES = '.D10.3/, 

* 6X, 'RFLIQ FOR THE HE ORIFICE = ',D10.3/, 

* 6X, 'RFGAS FOR THE HE ORIFICE = '.D10.3/, 

* 6X, 'RFLIQ FOR THE HE VALVE = '.D10.3/, 
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* 6X , ' RFGAS FOR THE HE VALVE = '.D10.3/, 

* 6X,' INJECTOR AND MANIFOLD CD = '.D10.3/, 

* 6X,' UPSTREAM SOURCE PRESSURE = ’,D12.5,’ PA'/, 

* 6X,' INITIAL BACKFLOW INTERFACE POSITION = ',1 PD 10. 3) 

1500 FORMAT(/4X, 'PIPE GEOMETRY:'/, 

* 6X, ' PIPE ',2X,' NODE l',2X,'NODE 2',2X,'AREA (M**2)', 

* 2X, 'LENGTH (M) ' / , (7X,I2, 5X, I2,6X,I2,5X, 1PD10.4, 

* 3X.0PF10.5)) 

1600 FORMAT ( /4X, ' PIPE INITIAL CONDITIONS: '/ ,6X, ' PIPE' ,2X, 

* 'VELOCITY (M/SEC )',2X,' LIQUID/ VAP INT. POSITION', 

* ’ (% LENGTH) ' / , ( 7X,I2,5X, 1PD11 .4,4X, 1PD11 .4) ) 

1700 F0RMAT(/4X, 'INITIAL NODAL DENSITIES: '/, 6X, 'NODE' ,2X, 

* 'DENSITY (KG/M**3) ' / , (7X,I2,6X, 1PD11 .4)) 

1800 F0RMAT(12(A4,2X)) 

1900 FORMAT (/2X, ****** SPECIES PROPERTIES *****’//, 

* 4X, 'MAJOR CONSTITUENTS :',/,6X,' INDEX', 2X,’ NAME', 

* 4X, 'MOLEC. WT. ' ,/,(8X,I2,3X,A4,4X,lPD11.4)) 

2000 FORMAT(/4X, 'RATE DETERMINING RADICALS',/, 

* 6X, 'INDEX', 2X, 'NAME', 4X, 'MOLEC. WT. ',/, 

* (8X,I2,3X,A4,4X,1PD11.4)) 

2100 FORMAT ( /2X, ' ***** "Z" ARRAY OF THERMO PROPERTY POLYNOMIAL ' 

* 'COEFFICIENTS *****',// ,6X,4HNAME,2X, 12HCOEFFICIENTS, // , 

* (6X,A4,5(2X, 1PD12.5)/ , 10X,5(2X,1PD12.5)/ ,10X, 

* 4(2X, 1PD12.5)//)) 

2200 FORMAT (/2X '***** REACTION RATE INDICIES AND DATA *****'//, 

* 4X, 'REACTION INDICIES ',//, 6X , 7HRXN ( J) ,4X,7HID( 1 , J) , 

* 2X, 7HID( 2 , J) , 2X, 7HID( 3 , J) , 2X, 7HID( 4 , J ) , / , 

* (9X.I2.8X, 4(12, 7X))) 

2300 FORMAT(/4X, 'ISIDE ARRAY' ,//,6X,3HRXN,4X,12HISIDE(K,RXN),/, 

* (7X,I2,4X,7(2X,I2))) 

2400 FORMAT(/4X, 'REACTION RATE DATA IN SI UNITS',//, 

* 6X , 3HRXN , 2X , 4HMODE , 1 5X , 2HBX , 1 OX , 3HTEN , 8X , 4HTACT , / ) 

2450 FORMAT( 7X,I2, 3X.I1 ,4X, 3HFWD, 2X,F11 . 3,2X,F11 .3,2X,F11 . 3/ , 

* 17X,3HREV,2X,F11 . 3,2X,F11 . 3 ,2X,F11 .3) 

2500 FORMAT (/2X , ' ***** CC AND INLET INITIAL CONDITIONS *****’//, 

* 4X, 'SPECIE', 2X, 'NAME', 2X,'S2 (KGMOLE/KG) ' ,2X, ' ISTRM’ , 2X, 

* 'SI (KGMOLE/KG)' ,2X,' STATE’,/, 

* (6X,I2,4X,A4,3X, 1PD12. 5,5X,I1 ,5X,1PD12.5,4X,A4)) 

2600 FORMAT(/4X, 'CC TEMPERATURE = '.0PF12.4,’ K'/, 

* 4X, 'CC PRESSURE = ',1PD12.4,' PA'/, 

* 4X, ' CC DENSITY = '.1PD12.5,' KG/M**3’/, 

* 4X, ' AVG. MOL. WT. = '.1PD12.5,' KG/KGMOLE' / , 

* 4X, ' CC VOLUME = '.1PD12.5,' M**3 ',//, 

* 4X , 6HSTREAM , X , 4HNAME , 2X , 4HTEMP , 8X , 4HMDOT , 9X , 2HSF , 1 2X , 2HHS , 

* 12X,2HSH/,18X,3H(K),6X,8H(KG/SEC),3X,12H(KGMOLE/SEC),3X, 

* 10H( J/KGMOLE) ,3X, 7H( J/SEC)/ / , (6X,I1 ,4X,A4, 2X.0PF6. 2,2X, 

* 1PD11 .5,2X, 1PD11 . 5, 2X,D12. 5,2X,D12.5)) 

2700 FORMAT(2X, 'TIME (MSEC) ' ,2X, 'MOLE NUMBER (KGMOLES/KG *1.0D+06) 

* ') 

2800 FORMAT(2X,F8. 3 , 7(X, 1PD9 . 2) ) 


2750 FORMAT (4X,' TIME \3X,' TEMPERATURE ' ,X, ' PRESSURE', X, ' OX MDOT 

* 2X, ' LIQ DEN' ,4X, 'DROP DIAM ',2X,'VAP RATE’ /, 3X, ' (MSEC) 6X, 

* '(K)’ ,6X,'(MPA)' ,4X,' (KG/S)' ,4X, ' (KG/M3) ' ,4X,' (MICRONS) ’ ,2X, 

* ' (KG/M3 S)') 

2850 FORMAT(2X,F8.3,2X,F8.2,X,F8.3,2X,F8.4,3X,lPD9.2,3X, 

* 0PF8.3.3X.F8.2) 


o o 


FILE INIT: INITIALIZATIONS OF VARIABLES 

IMPLICIT REAL* 8 (A-H.P-Z) 

IMPLICIT INTEGER* A (I-O) 

REAL*8 OMV 

CHARACTER*A ANAM(20) ,AST(20) , ASTNST( 10) ,ANST( 10) ,ALIQ,AGAS 


o o 


COMMONS S . PROC 


C 


COMMON 

* /CPARA/ SMW( 10) ,AST,ANAM,RATE( 10) , YPP( 10) ,QMV,FMV, PMV, 

* EMV , SMEMV , EGAS , RGAS IN , VOL , Q , S2MIN , S 2MLG , SM , SMINV , 

* TMSD , RGMX , ALIQ , AGAS , VR , RHOGAS , RHOLI Q , RH0H2 , CPSUM , 

* AH2IN, A02IN , VH2 , DM , DMMIN , DMC , CFTP , PHG , DPIN , RPX , 

* DLIQMN,FVI2MN,XPER,TMPER,PAMAX,CVR, 

* / CSVAR/ S2 ( 10 ) , TK , TKINV , PA, RHOCC , RHINV , DLIQ , XMASS , 

* /CSTRM/ SI (3, 30) ,ISTRM( 10) ,TKS( 3) ,EMS(3) ,DHF( 10) , 

* INST(10),HS(3),ASTNST,ANST,SH(3),SF(3), 

* /CINTE/ TMEND, TMI.TMPRNT, STEP, EPS I, METH, MITER, INDEX, 

* IER.NDGEAR, 

* / CIDX1 / NDEBUG(20) ,NREAC,NPROD,NRAD,NSRATE, JJ.NSTRM, 

* NV,NS,NS1,NS2,NS3,NS4,NS5,NS6,NS7,NS8,NS9,NS10,MCON, 

* KINET, 

* /CRATE/ ID(4,20) ,ISIDE( 10,20) ,MODR(20) ,BX(20) ,BX2(20) , 

* TEN(20) ,TEN2(20) ,TACT(20) ,TACT2( 20) ,ARR( 10) ,RIP( 10,10), 

* Xl(20) ,X2(20) , 

* / CPROP/ Z(200) ,H( 10) ,CP( 10),U(10) ,S( 10) ,DCP( 10), HO (10), 

* CP0(10),S0(10), 

* / CSSRR/ EPSISS,RLX,ITMAXS,IDGT, JER, 

* /CPIPE / V( 15) ,FVI( 15) ,RHO( 15) ,P( 15) ,AREA( 15) ,PL( 15) , 

* PS , FVI2 , CDIN , SM1MIN , RFLIQ , RFGAS , RFLV , RFGV , RFGHO , RFLHO , 

* / CIDX2/ Nl( 15) ,N2( 15) ,NN,NP,NHEO,NFVI,NPV 
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SUBROUTINE CSLOPE ( N , T , Y , YPRIME ) 


C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


CSLOPE CALCULATES THE DERIVATIVES WITH RESPECT TO TIME OF 
THE SYSTEM VARIABLES. 

N IS THE TOTAL NUMBER OF DEPENDENT VARIABLES. 


YPRIME( 1 :NS) 
YPRIME(NSl) 
YPRIME (NS2) 
YPRIME(NS3) 
YPRIME(NS4:NS5) 
YPRIME (NS 6: NS 7) 
YPRIME (NS 8: NS 9) 
YPRIME(NSIO) 


--> D(L0G(S2(I))/DT 
--> D(DLIQ)/DT 
--> D(LOG(TK))/DT 
--> D(LOG(PA))/DT 
--> D(V( J) )/DT 
--> D(FVI(J))/DT 
— > D(RHO(K))/DT 
— > D(FVI2)/DT 


SYSTEM VARIABLES INCLUDE COMBUSTION CHAMBER VARIABLES AND 
FEED SYSTEM VARIABLES: 


S2(I) IS MOLE NUMBER OF SPECIES I (KGMOLES I/KG TOTAL) 
DLIQ IS THE COMBUSTION CHAMBER LIQUID DENSITY (KG/M**3) 
TK IS COMBUSTION CHAMBER TEMPERATURE (K) 

PA IS COMBUSTION CHAMBER PRESSURE (PA) 

V(J) IS FLUID VELOCITY IN PIPE J. 

FVI(J) IS POSITION OF THE LIQUID/VAPOR INTERFACE IN PIPE J. 
RHO(K) IS DENSITY OF NODE K. 

FVI2 IS THE POSITION OF THE BACKFLOW INTERFACE (FROM THE 
COMBUSTION CHAMBER INTO THE OXIDIZER FEED SYSTEM. 


OTHER VARIABLES: 

OMV = VOLUMETRIC MASS FLOW RATE THROUGH THE OXIDIZER 
FEED SYSTEM 

FMV = VOLUMETRIC MASS FLOW RATE THROUGH THE FUEL FEED 
SYSTEM (CONSTANT) 

PMV = VOLUMETRIC MASS FLOW RATE OF HELIUM INTO THE 
COMBUSTION CHAMBER. 

EMV = EXIT VOLUMETRIC MASS FLOW RATE (THE THE EXIT 
TURBINE 

SM = AVERAGE MOLECULAR WEIGHT IN THE COMBUSTION CHAMBER 
(DOES NO TAKE LIQUID SPECIES INTO ACCOUNT) 

INCLUDE 'INIT. FOR/LIST ’ 

DIMENSION Y(N),YPRIME(N),ALG(20) 

\ 

INCLUDE 'COMMONSS. FOR/ LIST' 


--- RESET YPRIME ARRAY 
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c 

c 

c 


c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


DO 100 1=1 ,N 

100 YPRIME(I) =0.0 

--- IF THE HELIUM/ OXYGEN INTERFACE IS AT THE LAST NODE STOP 
EXECUTION 

IF(FVI(NP).GE.1.0) THEN 
WRITE (30,1000) 

CALL OUTPT(Y,TM,N,NSTEPS) 

STOP 

- ENDIF 

COUNT ITERATIONS AT EACH TIME STEP 

IF(T.EQ.TMOLD) THEN 
ITER = ITER + 1 
ELSE 

ITER = 1 
TMOLD = T 
ENDIF 

--- COMPUTE MOLE NUMBERS. NO SPECIES MOLE NUMBER MAY EXCEED 
1.0/SMW(I) (SMW(I) IS THE MOLECULAR WEIGHT OF SPECIES I) 

TMRAT =0.0 
SM = 0.0 
DO 110 1=1, NS 

ALG(I) = DLOG(1.0/SMW(I)) 

IF(Y(I) .LT.ALG(I) .AND.Y(I) .GT.S2MLG) THEN 
S2(I) = DEXP(Y(I)) 

ELSE IF(Y(I) .GE. ALG(I) ) THEN 
S2(I) = DEXP(ALG(I)) 

Y(I) = ALG(I) 

ELSE 

S2(I) = S2MIN 
Y(I) = S2MLG 
ENDIF 

TMRAT = TMRAT + S2(l)*SMW(I) 

110 CONTINUE 

— NORMALIZE MOLE NUMBERS 

DO 120 1=1, NS 

SM = SM + S2(I) 

S2(I) = S2(I)/TMRAT 

120 CONTINUE 

SMINV = l./SM 

--- FIND TEMPERATURE, PRESSURE, DENSITY AND THE LIQUID DENSITY 
FOR THE COMBUSTION CHAMBER AND BACKFLOW INTERFACE POSITION. 
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c 

c 


c 

c 

c 

c 

c 


DLIQ = Y(NS1) 

IF ( DLIQ . LT . DLIQMN ) THEN 
DLIQ = DLIQMN 
ENDIF 

TK = DEXP(Y(NS2) ) 

IF(TK.GT. 1 . 0D+04) THEN 
TK = 1.0D+04 

ELSE IF(TK.LT. 100.0) THEN 
TK = 90.0 

ENDIF ~~~ 

PA 5_DEXP(Y(NS3)) 

IF ( PA . GT . PAMAX ) THEN 
PA - PAMAX 

ELSE IF(PA.LT. 1 .0D+03) THEN 
PA = 1.0D+03 
ENDIF 

TKINV = 1.0 / TK 
FVI2 = Y(NS10) 

RHOCC = PA*RGASIN*SMINV*TKINV 
IF(RHOCC.LT.O.Ol) RHOCC = 0.01 
RHINV = 1.0 /RHOCC 

— - FILL ARRAYS FOR V, RHO AND FVI 

DO 130 1=1, NP 

V(I) = Y(NS3+I) 

FVI(I) = Y(NS5+I) 

130 CONTINUE 

DO 140 1=1, NN 

RHO(I) = Y(NS7+I) 

140 CONTINUE 

--- DETERMINE MASS FLOW RATE IN THE OXIDIZER FEED SYSTEM AND 
RATE OF ENERGY INFLUX TO THE CHAMBER FROM THE OXIDIZER 
FEED SYSTEM. 

IF(FVI(NP).LT.1.0) THEN 
PMV =0.0 
SF(3) = 0.0 
SH(3) = 0.0 

IF(V(NP) .GT.O.O.AND.RHO(NN) .GT.RGMX) THEN 
OMV = RHOLIQ * AREA(NP) * V(NP) / VOL 
SH(1) = OMV * HS( 1 ) / SMW(3) 

SF( 1) = 0.0 
ELSE 

OMV = RHOCC * AREA(NP) * V(NP) / VOL 

SH( 1) = 0.0 

SF( 1 ) = 0.0 

CALL THERPS(NS.TK) 
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DO 150 1=1, NS 

SH( 1) = SH( 1) + 0MV*S2(I)*H(I) 

150 CONTINUE 

ENDIF 
ELSE 

• OMV = 0.0 

PMV = RH0(NN)*AREA(NP)*V(NP)/V0L 
SF(3) = PMV/SMW( 1) 

SH(3) = SF(3)*HS( 1) 

ENDIF 

--- DETERMINE THE VOLUMETRIC EXIT MASS FLOW RATE (EMV) 

PR = PHG/PA 

PB = PR**1 .4286 - PR**1 . 7143 
IF(PB.LE.O.O) THEN 
EMV =0.0 
ELSE 

EMV = CFTP*DSQRT(PA*PB/TK) 

ENDIF 

EMV = EMV/ VOL 
IF (EMV. LT. 0.0) EMV = 0.0 

--- COMBRT CALCULATES THE RATE OF PRODUCTION OF EACH SPECIES IN 
THE CHAMBER DUE TO COMBUSTION AND VAPORIZATION. 

CALL COMBRT (T,Y,YPRIME) 

--- LOOPS TO CALCULATE YPRIMES C C --- YPP(I) ARE THE TERMS USED 
TO FIND D(TEMP)/D(TIME) 

DO 160 IB = 1,10 
160 YPP(IB) = 0.0 

--- CALL THERPS FOR THE SPECIFIC MOLAR INTERNAL ENERGY AND 
ENTHALPY FOR EACH SPECIES C 
CALL THERPS ( NS, TK) 

--- CONSERVATION OF MASS FOR THE GAS PHASE 

IF(V(NP) .GE.O.O.AND.RHO(NN) .GE.RGMX) THEN 
SUMM = VR + FMV + PMV 
ELSE 

SUMM = VR + OMV + FMV + PMV 
ENDIF 

RPX = SUMM - EMV 

--- OUTPUT VALUES TO THE SCREEN IF NDEBUG(IO) = 1 
IF(NDEBUG( 10) ) THEN 
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WRITE (6,1100) T,V(NP) ,TK,VR,DM,DLIQ 
ENDIF 
C 

C CONSERVATION OF SPECIES 

C 

• C NOTE: IF THERE IS BACKFLOW INTO THE PIPES OR THERE IS 

C AN INTERFACE IN THE PIPES DUE TO BACKFLOW (FVI2>0.0), 

C OMV ACCOUNTS FOR THE INFLOW OR OUTFLOW OF CHAMBER CONTENTS. 

C CHAMBER GAS IN THE PIPES IS ASSUMED TO HAVE THE SAME 

C DENSITY AND TEMPERATURE AS THE CHAMBER CONTENTS. 

C 

DO 170 1=1, NS 

IF(OMV.LT.O.O.AND.RHO(NN) .LT.RGMX) THEN 
IF(ISTRM(I).GT.O) THEN 

YP = (SF(ISTRM(I))-(SUMM-0MV)*S2(I)+RATE(I))*RHINV 
ELSE 

YP = (-(SUMM-0MV)*S2(I)+RATE(I))*RHINV 
ENDIF 
ELSE 

IF(ISTRM(I) .GT.O) THEN 

YP = (SF(ISTRM(I))-SUMM*S2(I)+RATE(I))*RHINV 
ELSE 

YP = ( -SUMM*S2(I)+RATE(I) )*RHINV 
ENDIF 
ENDIF 

YPRIME(I) = YP/S2(I) 

IF(Y(I) .GE.ALG(I) .AND. YP.GT. 0. 0) THEN 
YPRIME(I) = 0.0 
YP = 0.0 
ENDIF 

YPP(A) = YPP(4) + EMV*S2(I)*H(I) 

YPP(5) = YPP(5) + U(I)*S2(I) 

YPP(6) = YPP(6) + U(I)*YP 
YPP(8) = YPP(8) + S2(I)*(CP(I)-RGAS) 

YPP(9) = YPP(9) + YP 
170 CONTINUE 
YPP(l) = SH( 1 ) 

YPP(2) = SH( 2) 

YPP(3) = SH( 3) 

YPP(5) = YPP(5)*RPX 
YPP(6) = YPP(6)*RHOCC 
YPP(7) = Q/VOL 
YPP(8) = YPP(8)*RHOCC 

CONSERVATION OF ENERGY 

TP = (YPP( 1 )+YPP(2)+YPP( 3)-YPP(4)-YPP( 5)-YPP(6) -YPP( 7 ) )/YPP(8) 
YPP( 8) = YPP(8) * TK 
YPRIME(NS2) = TP/TK 
C 
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--- CONSERVATION OF MASS FOR THE LIQUID PHASE 

IF(ASTNST(1).EQ.AGAS) THEN 
YPRIME ( NS 1) = 0.0 
ELSE 

IF(V(NP) .GE.O.O.AND.RHO(NN) .GE.RGMX) THEN 
YPRIME(NSl) = OMV-VR-EMV*DLIQ/RHOCC 
ELSE 

YPRIME(NSl) = -VR-(EMV-OMV)*DLIQ/RHOCC 
ENDIF 
ENDIF 

THE IDEAL GAS LAW 

YPRIME(NS3) = YPRIME(NS2) + RPX * RHINV + YPP(9) * SMINV 

--- SUBROUTINE PIPES DETERMINES YPRIME FOR THE PIPING SYSTEM 
DEPENDENT VARIABLES. PIPES IS CALLED IF MCON = 1. OTHERWISE 
OXIDIZER FEED RATE REMAINS CONSTANT 

IF (MCON) CALL PIPES (Y, YPRIME, T) 

--- WRITE A SUMMARY OF THE TERMS USED TO CALCULATE DERIVATIVES 

IF(NDEBUG(8).EQ.l) THEN 

WRITE(30,1200) T,ITER,TK,PA,RHOCC,EMV,DLIQ,VR, 

* ( I , EMS (l),SF(I), SH( I), 1=1, 3) 

WRITE( 30,1300) (IS,H(IS),U(IS),CP(IS),RATE(IS),IS=1,NS) 
WRITE(30,1400) ((IB,YPP(IB)),IB=1,9) 

ENDIF 

--- WRITE A SUMMARY OF YPRIME VALUES CALCULATED IN CSLOPE 

IF(NDEBUG(9).EQ.l) THEN 
WRITE(30, 1500) T,ITER 
WRITE(30,1600) 

WRITE(30 , 1700) ((I,S2(I),Y(I),YPRIME(I)),I=1,NS), 

* NS1 ,DLIQ, Y(NS1) , YPRIME(NS1 ) ,NS2,TK,Y(NS2) , 

* YPRIME(NS2),NS3,PA,Y(NS3),YPRIME(NS3),NS10, 

* FVI 2, YPRIME (NS 10) 

WRITE(30,1800) ( J,V( J) ,YPRIME(NS3+J) ,FVI( J) , 

* YPRIME(NS5+J) , J=1 ,NP) 

WRITE(30, 1900) (K,RHO(K) , YPRIME (K+NS 7) ,K=1 ,NN) 

ENDIF 

RETURN 

1000 FORMAT(2X, ****** INTERFACE HAS REACHED THE LAST NODE *****'//, 

* 2X, ' ***** EXECUTION HALTED ***** ' / / ) 

1100 FORMAT(X, ' T ’ , 1PD8 . 2,X, ' V, TK, VR, DM, DLIQ *,D8.1,X, 




* 0PF6. 1 ,3(X, 1PD8 . 1) ) 

1200 FORMAT (/2X, ****** DEPENDENT VARIABLE DERIV. PARAMS *****'//, 

* 4X, 'TIME = ',1 PD 12. 5, 2X,' ITERATION = \I5/,4X, 

* 'TEMPERATURE = ' ,0PF15. 4, 2X, ' PRESSURE = ' , 1PD12. 5/ ,4X, 

* 'DENSITY = ' , 1PD12.5,2X, 'EXIT VOLUMETRIC MASS FLOW = ', 

* 1PD12.5/,4X,'LIQ DENSITY = '.1PD12.5,' (KG/M**3) 7 ,4X, 

* 'VAPORIZATION RATE = ',1PD11.4,' KG/S M**3'//,4X, 'STREAM' ,2X, 

* ' MASS FLOW ' ,2X, ' MOLE FLOW ', 2X, ' ENTHALPY FL0W’/,14X, 

* ' (KG/SEC) ' ,2X, ' (KGMOLE/SEC) ' , 3X, ’ (JOULES/SEC) 'll , 

* 3(7X,I1,4X, 1PD11 .4, 2X, 1PD11 . 4,2X, 1PD11 .4/ )) 

1300 F0RMAT(/4X, 'IN CSLOPE THERMO PROPS AND RXN RATES ARE:'/, 

* 4X, ' I' ,2X, ' H ' ,2X, ' U ' ,2X, ' CP 

* 4X, ' A ' , / , (4X,I2,4(2X, 1PD12.5))) 

1400 FORMAT(/4X, 'TEMPERATURE DERIVATIVE TERMS ARE: '/, 

* (4X, ' YPP( 'll,') = ' ,1PD12.5,2X, 'YPP( 'll, ' ) = '.1PD12.5)) 

1500 FORMAT(/ /2X, ' ***** CLSOPE YPRIME ESTIMATES *****'//, 

* 4X , ' TIME = ' ,1PD12. 5, 2X,' ITERATION = ',15,//) 

1600 FORMAT (4X, 3H I ,2X,12H S2(I) ,2X,12H Y(l) ,2X, 

* 13H YPRIME(I) ,//) 

1700 FORMAT(4X,I3,2X, 1PD12 . 5, 2X, 1PD12.5,2X, 1PD12.5) 

1800 FORMAT(/4X, 'PIPE' ,2X, ' VELOC. (M/S) ’ ,2X, 'D(VELOC. )/D(T) ’ , 

* 2X, 'FL/V INTERFACE', 2X,'D(FVI)/D(TIME)'//, 

* (5X,I2,X,4(5X,1PD11 .4))) 

1900 F0RMAT(/4X, 'NODE', 2X, 'DENSITY (M**3) ' , 2X, ' D(RHO)/D(TIME) ' // , 

* (5X,I2,X,2(5X,1PD11.4))) 

END 
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SUBROUTINE COMBRT ( TM , Y , YPRIME ) 


COMBRT FINDS THE RATE OF PRODUCTION OF EACH OF THE SPECIES IN 
THE COMBUSTION CHAMBER DUE TO REACTION RATE AND VAPORIZATION. 


C 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 

c 


INCLUDE 'INIT. FOR/LIST ' 

DIMENSION Y(NV),YPRIME(NV) 

INCLUDE ' COMMONSS . FOR/LIST ' 

— FIND THE VAPORIZATION RATE, VR (KG/M**3 S) AND YPRIME(NSl) 

CPSUM = THE AVERAGE SPECIFIC HEAT AT CONSTANT PRESSURE FOR 
THE COMBUSTION CHAMBER GASES. 

XK = THERMAL CONDUCTIVITY OF HYDROGEN (LEAST SQUARES LINEAR 
FIT FOR XK AS A FUNCTION ONLY OF TEMPERATURE) 

XMU = VISCOUSITY OF HYDROGEN (LEAST SQUARES LINEAR FIT FOR 
XMU AS A FUNCTION ONLY OF TEMPERATURE) 

DATA PI/ 3. 141592654/ 

TAV = (TKS(l) + TK)/2.0 
CALL THERPS(NS,TAV) 

CPSUM = 0.0 
DO 100 1=1, NS 

CPSUM = CPSUM + CP(I)*S2(I) 

100 CONTINUE 

CALL THERPS(NS,TKS(1)) 

XK = 5.282D-02 + 4. 1>17D-04*TAV 
XMU = 5.0767D-06 + 1 .4884D-08*TAV 

(B) FIND THE AVBERAGE DROPLET DIAMETER, DM (M), THE RATE OF 
VAPORIZATION, VR (KG/ M3 S), AND THE TOTAL NUMBER OF 
DROPLETS IN THE COMBUSTION CHAMBER, XNDR. 

THE EXPRESSION USED TO CALCULATE EVAPORATION RATE WAS 
ADAPTED FROM WEBBER (1972) (J. OF THE ARS) 

IF(ASTNST(1) .EQ.AGAS) THEN 
VR = OMV 
GO TO 110 
ENDIF 

DM = DMC*OMV / FMV*DSQRT ( SMW ( 2 ) *AH2IN*PA/ RGAS / RHOLIQ/ TKS ( 2 ) ) 
IF(DM.LT.DMMIN) DM = DMMIN 
DM3 = DM** 3 

REN = RHOCC*VH2*DM/XMU 

ALPH = 1.0 + CVR * (TK-TKS(l)) 

IF( ALPH.GT. 1 . 0) THEN 

VRD = PI*XK*DM*( 2 . 0+0 . 5*REN**0 . 5 )*DLOG( ALPH) / CPSUM 
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VRL = 0.0 
ENDIF 

XNDR = 6 . 0*DLIQ*VOL/RHOLIQ/DM**3/PI 
VR = 6 . 0*VRD*DLIQ/PI/RHOLIQ/DM3 
IF(VR.LE.O.O) THEN 
VR = 0.0 
ENDIF 

110 CONTINUE 
C 

C DETERMINE THE RATES OF PRODUCTION DUE TO COMBUSTION AND 

C VAPORIZATION 

C (A) IF KINET = 1, CHEMICAL KINETICS ARE CONSIDERED. REACTION 
C RATE IS CALCULATED IN RADSS ASSUMING THAT ATOM AND RADI- 

C CAL CONCENTRATIONS ARE STEADY STATE. 

C (B) IF KINET = 0, ALL OXYGEN WHICH VAPORIZES IS ASSUMED TO 
C IMMEDIATELY COMBUST. 

C 

IF(KINET.EQ.l) THEN 
CALL RADSS 

RATE(3) = RATE(3) + VR/SMW(3) 

ELSE IF(KINET.EQ.O) THEN 
RATE(3) = 0.0 
RATE(4) = 2.0*VR/SMW(3) 

RATE( 2 ) * -2.0*VR/SMW(3) 

ELSE IF(KINET.EQ.2) THEN 
CALL RATES 
DO 200 1=1, NS 

RATE(l) = ARR(I) 

200 CONTINUE 

RATE(3) = RATE(3) + VR/SMW(3) 

ENDIF 

C 

IF(NDEBUG( 15) ) WRITE( 30 , 1000) TM , TK , PA , DM , VRD , XNDR , REN , VR , 

* (ANAM(I),RATE(I),I=1,NS) 

C 

RETURN 

C 

1000 FORMAT(/,2X, ****** OUTPUT FROM COMBRT ****** //,4X, ’ AT TIME 

* 1PD11.5/, 

* 4X, 'TEMPERATURE = ',0PF10.3,' K 7, 

* 4X, ' PRESSURE = ',1PD12.5,' PA'/, 

* 4X, 'DROPLET DIAMETER = ’,D12.5,' M '/, 

* 4X, 'DROPLET EVAPOATION RATE = ',D12.5,' KG/S/DROPLET'/, 

* 4X, 'NO. OF DROPLETS = ',D12.5,' DROPLETS'/ 

* 4X, 'REYNOLDS NO. OF DROPS = \D12.5/, 

* 4X, 'TOTAL EVAPORATION RATE = ’,D12.5,’ KG/M**3 S'//, 

* 4X, 7HSPECIES,2X, 'RATE OF PRODUCTION (KGMOLES/M**3 S) ' // , 

* (5X,A4,5X, 1PD12.5)) 


c 


END 
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SUBROUTINE DTSTR(NV , CSLOPE , TM , STEP , Y , TMNXT , WK , NDEBUG ) 


C 

C THIS IS A RUNGE-KUTTA FOURTH ORDER INTEGRATING SUBROUTINE 
C WHICH ANSWERS A CALL MIMICKING A CALL TO IMSL SUBROUTINE 
C DGEAR. THE PURPOSE OF DGEAERX IS TO ALLOW DEBUGGING OF 

C THE TRANSIENT CHUGGING CODE AT MSFC, WHERE DGEAR IS NOT 

C AVAILABLE. 

C 

C THE SIGNIFICANT VARIABLES ARE 

C NV=NUMBER OF VARIABLES 

C CSLOPE=NAME OF ROUTINE TO EVALUATE 1ST DERIVATIVE 

C TM=INDEPENDENT VARIABLE 

C STEP=DELTA TM 

C Y=ARRAY OF INDEPENDENT VARIABLES 

C SHOULD BE AN INTEGER 

C TMNXT=ENDING TM -( TMNXT -TM)/ STEP 

C WK=WORK VECTOR OF 4*NV SIZE 

C OTHERS ARE UNUSED DUMMY PARAMETERS 

C 
C 

IMPLICIT REAL* 8 (A-H.P-Z) 

DIMENSION Y(NV),WK(4,NV),YEND(100),F(100),NDEBUG(20) 

C 

C 

NSTEPS=( TMNXT -TM)/ STEP 
IF(NSTEPS.LE.O) NSTEPS=1 
C 
C 

DO 100 J=1,NSTEPS 
C 

C GET FIRST ESTIMATE OF THE DELTA X'S 
C 

CALL CSLOPE (NV,TM,Y,F) 

DO 10 1=1, NV 

WK(1,I)=STEP*F(I) 

YEND(I)=Y(I)+WK(1,I)*0.5 
10 CONTINUE 

C 

C GET SECOND ESTIMATE (USING YEND) 

C 

NE=2 

TM2E=TM+STEP*0 . 5 

CALL CSLOPE ( NV , TM2E , YEND , F ) 

DO 20 1=1, NV 

WK( 2 , I )=STEP*F ( I ) 

YEND(I)=Y(I)+WK(2,I)*0.5 
20 CONTINUE 

C 

C THIRD ESTIMATE 

C 
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NE=3 

CALL CSLOPE ( NV , TM2E , YEND , F ) 

DO 30 1=1, NV 

WK( 3 , I )=STEP*F ( I ) 

YEND (I)=Y(I)+WK(3,I) 

30 CONTINUE 

C 

C LAST ESTIMATE 
C 

NE=4 

TM2E=TM+STEP 

CALL CSLOPE(NV , TM2E , YEND , F ) 

DO 40 1=1, NV 

WK( 4 , I )=STEP*F ( I ) 

40 CONTINUE 

C 

C COMPUTE Y'S AT TM + STEP BY WEIGHTED AVERAGE 
C 

DO 50 1=1, NV 

YEND(I)=Y(l)+(WK(l,I)+2*WK(2,I)+2*WK(3,I)+WK(4,I))/6.0 
50 CONTINUE 

C 

C WRITE VALUES OF THE DERIVATIVE ESTIMATES AND "OLD" AND "NEW" 
C VALUE OF Y IF NDEBUG( 15) = 1 
C 

IF(NDEBUG( 15) ) THEN 

WRITE(30,1000) TM,STEP 

WRITE(30, 1100) (K, (WK(I,K) ,1=1 ,4) ,Y(K),YEND(K) ,K=1 ,NV) 
ENDIF 
C 

C ADVANCE Y&TM FOR NEXT STEP 
C 

DO 60 1=1, NV 
Y(I)=YEND(I) 

60 CONTINUE 

TM=TM+STEP 

C 

100 CONTINUE 
C 

RETURN 

C 

C FORMATS 
C 

1000- FORMAT (//,2X, ****** DTSTR OUTPUT *****'//, 

* 6X, 'TIME = ' , 1PD12.5, 1 SEC’ ,2X, 'STEP = ',D12.5,' SEC'//, 

* 4X, ' INDEX \X,' ESTIMATE 1 ' ,X, ' ESTIMATE 2' ,X, 'ESTIMATE 3', 

* X,' ESTIMATE 4', X, ' Y OLD \X,' Y NEW '//) 

1100 FORMAT(6X,I2,2X,1PD10.3,X,D10.3,X,D10.3,X,D10.3,X,D10.3, 

* X.D10.3) 

C 
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END 
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SUBROUTINE NPT 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 


c 

c 

c 

c 


--- DEEINITION OF INPUT VARIABLES. 


ANAM(I) 

S2(I) 

XIN(I) 

XCC(I) 

AST(I) 

ISTRM(I) 

INST(J) 

ANST(J) 

ASTNST(J) 

DHF(J) 

HS( J) 
TKS(J) 

SF( J) 

SH( J) 


= HOLERITH NAME OF SPECIES I 
= MOLE NUMBER OF SPECIES I (KGMOLES OF 1/ 

KG (TOTAL)) 

= INLET MOLE FRACTION OF SPECIES I 
= INITIAL MOLE FRACTION OF SPECIES I. 

= PHASE OF SPECIES I 

= NUMBER OF THE INLET STREAM CONTAINING SPECIES I 
= NUMBER OF THE SPECIES CONTAINED IN REACTANT 
STREAM J 

= HOLERITH NAME OF THE SPECIES IN STREAM J 
= PHASE OF REACTANT STREAM J 
= LATENT HEAT OF VAP. OF LIQUID IN STREAM J 
= ENTHALPY OF THE INCOMING FLUID IN STREAM J 
= TEMPERATURE OF STREAM J (ASSUMED CONSTANT) 

= MOLE FLOW (KGMOLE/M'3 SEC) INTO THE CC 
= ENTHALPY FLOW (J/M‘3 SEC) INTO THE CC 


INCLUDE 'INIT. FOR/ LIST* 


CHARACTER*4 AXP( 10) ,AXR(10) ,ASTP( 10) ,DATA( 12) ,ASTR( 10) , 

* AT(7),CMNT,PHAZ, DTI, DT2,ATYP,AX,ASTX, BLANK, 

* AREAC , APROD , ARADI , ACOMM , AREVE , AGLOB , AHETE , 

* APYRO.ACGS, THIRD 

DIMENSION SMWR( 10) ,XIN(20) ,XCC(20) ,SMWP( 10) ,XCCP( 10) , 

* TX(50) ,TY(50) 


INCLUDE 'COMMONSS. FOR/ LIST’ 


DATA BLANK/' '/ , AREAC/ 'REAC' /, APROD/ ' PROD' /, ARADI/ 'RADI '/ , 

* S2/ 10*0 .0/ ,Sl/90*0 . 0/ ,XIN/20*0 . 0/ ,XCC/20*0 . 0/ , 

* ISTRM/10*0/ .NREAC/O/ .NPROD/O/ ,NRADI/0/ ,NS/0/ ,NSRATE/ 0/ , 

* ALIQ/'L ' / , AGAS / ’ G 7, LSI/ 1 / .ACOMM/ ' COMM'/ , 

* AREVE/ ' REVE ' / .AGLOB/ ’ GLOB ' / .AHETE/ ' HETE'/ , APYRO/ ' PYRO'/ , 

* ACGS/'CGS ' / .THIRD/ *M 7, XMAX/0. 00125/, XMIN/0. 0005/, 

* TENLN/ 2. 302585093/ 

--- READ RELEVANT SPECIES AND ASSIGN THEM AS REACTANT, PRODUCT OR 
RATE DETERMINING RADICAL. 


100 READ(20, 1000) ATYP,AX,SMWX,ISTRMX,XINX,XCCX,ASTX,DHFX,HNSTRM 
IF (ATYP.EQ. BLANK) GO TO 110 
IF (ATYP.EQ. AREAC) THEN 
NREAC = NREAC + 1 
ANAM(NREAC) = AX 
ANST(ISTRMX) = AX 
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ft 


ft 


t> 


m 


SMW( NREAC) = SMWX 
ISTRM(NREAC) = ISTRMX 
XIN( NREAC) = XINX 
XCC( NREAC) = XCCX 
AST (NREAC) = ASTX 
ASTNST ( I STRMX ) = ASTX 
DHF( ISTRMX) = DHFX 
HS( ISTRMX) = HNSTRM 
INST (ISTRMX) = NREAC 
ELSE IF ( ATYP . EQ . APROD ) THEN 
NPROD = NPROD + 1 
AXP(NPROD) = AX 
SMWP( NPROD) = SMWX 
XCCP(NPROD) = XCCX 
ASTP( NPROD) = ASTX 
ELSE 

NRAD = NRAD + 1 
AXR(NRAD) = AX 
SMWR(NRAD) = SMWX 
ASTR(NRAD) = ASTX 
END IF 
GO TO 100 
110 CONTINUE 
C 

C --- PUT ALL SPECIES INTO THE S 2 AND OTHER ASSOCIATED ARRAYS. 
C ORDER IS: FIRST REACTANTS, SECOND PRODUCTS AND THIRD 

C RADICALS. 

C 

NS = NREAC + NPROD 
NSRATE = NS + NRAD 
NS1 = NS + 1 
NS2 - NS1 + 1 

NS3 = NS1 + 2 

DO 120 1=1, NPROD 
II = NREAC + I 
ANAM(II) = AXP(I) 

SMW(II) = SMWP(I) 

XCC(II) = XCCP(I) 

AST (II) = ASTP(I) 

120 CONTINUE 

DO 130 1=1, NRAD 
II = NS + I 
ANAM(II) = AXR(I) 

SMW(II) = SMWR(I) 

AST (II) = ASTR(I) 

130 CONTINUE 

DO 140 1=1, NSRATE 

S2( I) = XCC(I)/SMW(I) 

IF(S2(l) .LT.S2MIN) S2(I) = S2MIN 
140 CONTINUE 
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DO 150 1=1 ,NREAC 

150 S1(ISTRM(I) ,1) = XIN(I)/SMW(I) 

DEFINE THE ONE DIMENSIONAL Z ARRAY WHICH HOLDS THE COEFFI- 
CIENTS OF THE POLYNOMIALS THAT GIVE THERMODYNAMIC PROPERTIES. 

160 READ (20,1100) (DATA(I),I=1,3) 

DO 170 1=1 .NSRATE 

IF(DATA(1) .EQ.ANAM(I)) THEN 

READ (20,1200) (Z( J) , J=(I-1)*14+1 ,1*14) 

GO TO 160 "" 

END IF 
170 CONTINUE 

IF (DATA(l).EQ. BLANK) THEN 
READ (20,1350) 

GO TO 180 
ENDIF 

READ (20,1300) 

GO TO 160 
180 CONTINUE 


- SET INITIAL MOLE AND ENERGY FLOWS (SF AND SH) 


* 


DO 190 J-1,3 

SF(J) - EMS( J)/VOL/SMW(INST( J) ) 

SH( J) = SF( J) * HS( J) 

190 CONTINUE 
C 

C --- READ AND DETERMINE CONSTANTS USED IN DETERMINING 
C CHEMICAL REACTION RATES. 

C 

C --- READ MECHANISM/ RATE DATA CARDS 

C COLUMNS 13 THROUGH 48 OF THE MECHANISM CARD ARE USED 

C AS A LABEL TO IDENTIFY THE KINETIC MECHANISM USED ON 

C THE SOLUTION OUTPUT. 

C 

C THE VARIABLE DTI (COLUMNS 73/76) IS USED AS A FLAG: 

C 

C CGS ---> 

C 
C 

C COMM ---> 

C REVE ---> 

C 

C GLOB — -> 

C PYRO ---> 

C HETE ---> 

C 

C THE HETEROGENEOUS REACTIONS ARE PROGRAMED IN CHAPTERS 
C 1 & 2 OF SUBROUTINE CORECT. BX(J), TEN(J), AND TACT(J) 


CGS UNITS, RATE CONSTANTS IN GM-MOLES, CM, 
SEC AND EACT IN (KCAL/GM-MOLE) OTHERWISE 
THE SI UNITS ( KG -MOLES, M, SEC) ARE ASSUMED. 
COMMENT CARD, FIRST 48 CHARACTERS PRINTED 
REVERSE RATE DATA, IN SAME UNITS AS FORWARD 
DATA 

GLOBAL RATE EXPRESSION DATA IN SI UNITS 
HETEROGENEOUS PYROLYSIS REACTIONS 
HETEROGENEOUS CHAR REACTIONS 
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C ARE READ AS BX(J), ALPHA( 1 OR 2), AND TACT(J) FOR PYROLY- 
C SIS AND BX( J) , TEN(J) AND TACT(J) FOR THE CHAR REACTIONS. 
C 

C DTI AMD DT2 (COL 73/80) CAN HAVE ANYTHING (COMMENTS) IF 
C ABOVE SIX WORDS ARE NOT REQUIRED. 

C TACT IS ACTIVATION TEMPERATURE, = EACT/ GASCON, DEG K 
C 

IF(NDEBUG(3)) WRITE (30,1500) 

DO 200 1=1,20 
TEN2(I)=0. 

BX2(I)=0. 

200 CONTINUE 
C 

JJ=1 

NPYROL=0 

NGL0B=0 

NHETER=0 

C 

210 READ (20,1600) (DATA(I) ,1=1 , 12) ,BX( JJ) ,TEN( JJ) , 

* TACT( JJ) ,DT1 ,DT2 

IF (DATA(l).EQ. BLANK. AND. DTI. NE.ACOMM) GO TO 400 

--- CHECK FOR COMMENT CARD 

IF (DTI. NE.ACOMM) GO TO 220 

IF (NDEBUG(3)) WRITE (30,1700) (DATA(I),I=1,12) 

GO TO 210 

--- SCAN TO INSURE ALL REACTANTS IN SPECIES LIST 

220 DO 240 IR=1 , 6 
I=2*IR-1 

IF (DATA(I).EQ. BLANK. OR. DATA(I).EQ. THIRD) GO TO 240 
DO 230 IS=1 .NSRATE 

IF (DATA(I) .EQ.ANAM(IS)) GO TO 240 
230 CONTINUE 

IF ( NDEBUG( 3 ) . NE . 1 ) GO TO 210 
WRITE (30,1800) 

WRITE (30,1900) (DATA(IDD),IDD=1,12) 

WRITE (30,2000) DATA(I) 

GO TO 210 
240 CONTINUE 
C 

C CHECK FOR REVERSE RATE DATA. NOTE: 

C 

C (1) ORDER OF CARDS MUST BE CORRECT 

C (2) UNITS OF REVERSE DATA ASSUMED SAME AS FORWARD DATA 

C 

IF (DT1.NE.AREVE) GO TO 250 
J = JJ-1 
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BX2( J) = BX(JJ) 

TEN2( J) = TEN(JJ) 

TACT2 ( J ) = TACT(JJ) 

IF (NDEBUG(3)) WRITE (30,2100) BX2 ( J ) , TEN2 ( J ) , 
* TACT2( J) ,DT1 ,DT2 

CONVERT BX2 FOR INTERNAL CALCULATIONS 

BX2( J) = BX2( J)*TENLN 

IF (LSI) GO TO 210 

BX2( J) - BX2( J)-TENLN*3 . 0 

TACT2( J) = TACT2( J)*1000 .0/ 1 . 987 

IF (MODR( J) .EQ.2) BX2(J) = BX2( J)-TENLN*3.0 

GO TO 210 

--- CHECK FOR UNITS 


250 LSI = 1 

IF (DT1.EQ.ACGS) LSI = 0 

IF (LSI.EQ.0) TACT( JJ)=TACT( JJ)*1000 . / 1 .987 
IF (NDEBUG(3) ) WRITE (30,2200) JJ, (DATA(I) ,1=1 , 12) , 

* DTI , DT2 , BX( J J) , TEN( JJ ) , TACT( J J ) 

IF (DTI .NE.AFYRO) GO TO 260 
NPYROL = NPYROL+1 
IF (NPYROL. GT. 2) WRITE (30,2300) 

NHETER = NHETER+1 
NGLOB = NGLOB+1 
260 CONTINUE 

IF (DT1.NE.AHETE) GO TO 270 
NGLOB = NGLOB+1 
NHETER = NHETER+1 
270 CONTINUE 

IF (DT1.EQ.AGLOB) NGLOB=NGLOB+ 1 

CONVERT BX FOR INTERNAL CALCULATIONS 

BX(JJ) = BX( JJ)*TENLN 

ID(I, J) IS THE INDEX NUMBER OF THE I-TH DISTINCT SPECIES 

IN REACTION J ... 1=1,4 AS NO DISTINCT THIRD BODIES ARE 
CONSIDERED 

9 > DO 280 1=1,4 

280 ID(I, JJ) = 0 
C 

ND = 1 

DO 330 N-1,6 
K = N*2 - 1 

^ IF (DATA(K).EQ. BLANK) GO TO 330 

IF (DATA(K).NE. THIRD) GO TO 290 
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DATA(K) = BLANK 
GO TO 330 
290 CONTINUE 

DO 300 1=1 .NSRATE 

IF (DATA(K) .NE.ANAM(I)) GO TO 300 
II = I 
GO TO 310 
300 CONTINUE 

310 IF (K.GT.3) GO TO 320 

-- ID(ND,JJ) = II 

ND = ND+1 
GO TO 330 

320 IF (ND.EQ.2) ND = 3 

ID(ND, JJ) = II 
ND = ND+1 
330 CONTINUE 
C 
C 

C — - STORE THE TYPE OF REACTION ... THREE TYPES 
C MODR 1 ... A + B — > C + D 

C MODR 2 ... AB + M ---> A + B + M 

C MODR 3 ... A + B + M ---> AB + M 

C 

MODR( JJ)=1 

IF (ID(2,JJ).EQ.0) MODR( JJ)=2 
IF (ID(4,JJ).EQ.0) MODR(JJ)=3 
C 

C — THE FOLLOWING SECTION, UP TO STATEMENT 355 INCLUSIVE, 

C MAY BE ELIMINATED IF REVERSE (AS WELL AS FORWARD) RATE 

C DATA IS SUPPLIED FOR ** ALL ** REACTIONS. 

C 

C THIS SECTION CALCULATES REVERSE RATE CONSTANTS FROM 
C EQUILIBRIUM CONSTANTS AND FORWARD RATE CONSTANTS FOR 
C TWENTY POINTS OVER THE TEMPERATURE RANGE 1000K TO 4000K 
C 

IF (DT1.EQ.AGLOB) GO TO 390 

IF (DT1.EQ.AHETE) GO TO 390 

IF (DT1.EQ.APYRO) GO TO 390 

DX=(XMAX-XMIN)/19.0 

SUMX=0 . 0 

SUMY=0 . 0 

IHCPS=2 

DO 360 N-1,20 

TX(N)=XMIN + DX*DFL0AT(N-1) 

SUMX=SUMX+TX(N) 

XTKINV=TX(N) 

XTK=1 . 0/XTKINV 
TLN=DLOG(XTK) 

CALL THERPS( NSRATE, XTK) 

RTKINV = 1 . / (XTK * RGAS) 
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SUM 1=0 . 0 
DO 350 ND=1,4 
K=ID(ND, JJ ) 

IF (K.EQ.O) GO TO 350 
GF = H0(K) - S0(K) 

IF (ND.LT.3) SUM1=SUM1+GF 
IF (ND.GE.3) SUM1=SUM1-GF 
350 CONTINUE 
TM1=0 . 0 

IF (ID(2,JJ).EQ.O) TMl=TLN-2. 50034 
IF (ID(4,JJ).EQ.0) TM1=2. 50034 -TLN 

-2.50034=LN(RGAS) RGAS IN UNITS ATM M**3/KGMOLE K HERE ONLY 

TY(N)=TM1-SUM1+TEN(JJ)*TLN-TACT(JJ)*XTKINV+BX(JJ) 

SUMY = SUMY+TY(N) 

360 CONTINUE 

XBAR = SUMX/20.0 
YBAR - SUMY/20.0 
SUMX = 0.0 
SUM1 = 0.0 
SUMY = 0.0 
DO 370 N=l,20 

SUMX = SUMX+TY(N)*(TX(N)-XBAR) 

SUM1 = SUM1+(TX(N)-XBAR)**2 
SUMY = SUMY+(TY(N)-YBAR)**2 
370 CONTINUE 

TEN2( JJ) = 0.0 

TACT2( JJ) = -SUMX/SUM1 

BX2( JJ) = ( YBAR+TACT2 ( J J ) *XBAR ) /TENLN 

SUMX = 0.0 

DO 380 N=l,20 

SUMX=SUMX+( TY ( N ) +TACT2 ( J J ) *TX( N ) -TENLN*BX2 ( J J ) ) **2 
380 CONTINUE 

SUMY = DSQRT( 1 . 0-SUMX/ SUMY) 

SUMX = DSQRT( SUMX/ 19.0) 

DATAT = TACT2 ( JJ ) 

IF (LSI.NE.l) DATAT=TACT2(JJ)*1. 987*0. 001 
IF (NDEBUG(3)) WRITE (30,2400) SUMX, SUMY, BX2( JJ) , 

* TEN2( JJ),TACT2( JJ) 

CONVERT BX2 FOR INTERNAL CALCULATIONS 

BX2( JJ) = BX2 ( J J ) *TENLN 

390 JJ = JJ+1 

--- CONVERT ALL RATE DATA TO SI UNITS 




IF (LSI) GO TO 210 
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J = JJ-1 

BX( J) =BX( J ) -TENLN*3 . 0 
BX2 ( J )=BX2 ( J ) -TENLN*3 . 0 
IF (MODR( J) .EQ.2) BX2( J)=BX2( J)-TENLN*3 . 0 
IF (MODR( J) . EQ. 3) BX(J)=BX( J)-TENLN*3.0 
GO TO 210 

400 JJ=JJ-1 

NPYROP=NPYROL+ 1 
NGL0BP=NGL0B+ 1 
NHETRP=NHETER+ 1 

--- SET CONTACT INDEXES TO UNITY 

DO 410 J=1,JJ 
X1(J)=1.0 
X2( J)=l .0 
410 CONTINUE 

RETURN TO TSTR 

RETURN 

1000 FORMAT(A4,X,A4,X,F7 .4,X,I1,X, 1PD11 . 5,X,D11 .5,X,A4,X, 

* D11.4,X,D11.4) 

1100 FORMAT ( 3 A4) 

1200 FORMAT(5D15.8/,5D15.8/,4D15.8/) 

1300 FORMAT( /// ) 

1350 FORMAT(/) 

1400 FORMAT( (A4,5(X, 1PD11 .4)) )' 

1500 FORMAT ( /2X, 1 ***** REACTION RATE DATA FROM NPT *****'//) 

1600 FORMAT (12A4.3F8. 3, 2A4) 

1700 FORMAT(6X,3H***,12A4,3H***) 

1800 FORMAT( / 4X, 26HF0LL0WING REACTION IGNORED) 

1900 FORMAT (8X, 6 A4,5H >,6X,6A4) 

2000 FORMAT (4X.9HREACTANT ,A4,26H NOT FOUND IN SPECIES LIST) 

2100 FORMAT(/14X, 17HREVERSE RATE DATA, 28X, 3F15 . 3,2A4) 

2200 FORMAT ( / 4X , 9HREACTI ON ,I2,3X,6A4,5H >,3X,6A4 /, 

* 6X, 14HDATA COMMENT: , 2A4/ , 6X.21HFWD RATE DATA: BX = , 

* F10.3,X,7HTEN = ,F10.3,X,8HTACT = ,F10.3) 

2300 FORMAT(2X, ****** TOO MANY PYROLYSIS REACTIONS ***** • ) 

2400 FORMAT ( 6X , 3 3HCALCULATED REVERSE RATE CONSTANTS/, 

* 6X, 16HSTD DEVIATION ■ , 1PD11 .3,X,13HCORR COEFF = ,D11.3/, 

* 6X.21HREV RATE DATA : BX2 = ,0PF10.3,X,7HTEN2 = ,F10.3,X, 

* 8HTACT2 = ,F10.3) 

2500 FORMAT(2X, ****** KINETIC REACTION RATE DATA IN SI UNITS ' 

* • ***** » // 1 

* 7X, 1HJ,6X, 4HM0DR , 1 2X , 2HID , 1 9X , 2HBX , 1 OX , 3HTEN , 9X , 4HTACT , 1 3X , 

* 3HBX2 , 9X , 4HTEN2 , 9X , 5HTACT2/ ) 

2600 FORMAT( (5X,I2, 1H. , 18, 3X,4I5,2( 3X, 3F13 . 3) )) 

END 


SUBROUTINE PIPES (Y.YPRIME.TM) 

C 

C 

INCLUDE 'INIT. FOR/ LIST* 

C 

DIMENSION Y(NV) ,YPRIME(NV) ,DV( 15) ,DRH( 15) ,FVIC( 15) , 

* SUM1( 15) ,SUM2( 15) ,RF( 15),ALEQ(50,50),BLEQ(50), 

* WKAREA(IOO) 

C 

• INCLUDE 'COMMONSS. FOR/LIST' 

C 

C SHUTDOWN IS BEGUN BY SETTING THE DENSITY OF THE FURTHEST 

C UPSTREAM NODE IN THE PIPING SYSTEM EQUAL TO THE DENSITY 

C OF THE HELIUM AT THE HELIUM SOURCE PRESSURE AND TEMPERA - 

C TORE. THIS IS DONE AT TIME = TMSD (TIME OF SHUTDOWN). 

C 

IF(TM.GT.TMSD) THEN 
RHO(l) = RHOGAS 
Y(NS8) = RHOGAS 
ENDIF 

IF(NDEBUG( 12) ) THEN 
WRITE(30, 1000) 

ENDIF 

C 

C — FIND THE LIQUID/ VAPOR INTERFACE POSITION AND VELOCITY IN 
C EACH PIPE. NOTE: 

C 

C (1) POSITION IS FIXED USING POSITION AND NODAL DENSITIES 
C (0.0 < FVI < 1.0). 

C (2) IF A NODE HAS JUST BECOME GASEOUS (THE INTERFACE HAS 
C JUST PASSED THROUGH) ITS DENSITY IS RECALCULATED USING 

C THE PREVIOUS NODAL PRESSURE AND THE IDEAL GAS RELATION. 

C (3) D(FVI( J))/DT IS EITHER ZERO (IF THERE IS NO INTERFACE 
C IN THE PIPE) OF EQUAL TO THE FLUID VELOCITY IN THE PIPE. 

C (4) AN INTERFACE CAN ENTER THE PIPE VIA BACKFLOW FROM THE 
C COMBUSTION CHAMBER. THIS OCCURS WHEN THE VELOCITY IN THE 

C LAST PIPE IS NEGATIVE. THE INTERFACE THAT RESULTS FROM 

C BACKFLOW IS NEVER ALLOWED TO PROCEED PAST THE UPSTREAM 

C NODE OF THE LAST PIPE. 

C 

DO 100 J=1 ,NP 
ND1 = Nl( J) 

ND2 = N2( J) 

IF(J.LE.NHEO) THEN 
RFLX = RFLHO 
RFGX = RFGHO 
ELSE IF(J.EQ.NPV) THEN 
RFLX = RFLV 
RFGX = RFGV 
ELSE 


171 


on non oonnonnnooooonn nnn 


RFLX = RFLIQ 
RFGX = RFGAS 
ENDIF 

BOTH PIPE NODES ARE GASEOUS 

IF(RH0(ND1) .LT.RGMX. AND. RH0(ND2) .LT.RGMX) THEN 
YPRIME ( NS 5+ J ) =0.0 
FVI(J) = 1.0 
FVIC(J) = 0.0 

RFJ = RFGX*PL(J)*V(J)*V(J)*RH0(ND2) 

RF( J) = DSIGN(RFJ,V(J>) 

Y(NS5+J) =1.0 

— ONLY THE UPSTREAM NODE IS GASEOUS. THREE CASES: 

(1) 0 < FVI < 1.0: THERE IS AN INTERFACE IN THE PIPE 

(2) FVI > 1.0: THE INTERFACE HAS JUST LEFT THE PIPE. 
THE DOWNSTREAM NODE IS MADE GASEOUS AND FVI IS SET 
AT 1.0. 

(3) FVI < 1.0: THE INTERFACE HAS RETURNED TO THE UP- 
STREAM PIPE. DENSITY OF NODE ND1 RETURNS TO THE 
CONSTANT LIQUID DENSITY (RHOLIQ) . FVI IS RESET TO 
0.0 AND FVI AND ITS VELOCITY ARE RECALCULATED FOR 
THE UPSTREAM PIPE. 

ELSE IF(RHO(NDl). LT.RGMX) THEN 

THE INTERFACE HAS JUST LEFT THE PIPE 

IF(FVI(J).GE.1.0) THEN 
YPRIME ( NS 5+ J ) =0.0 
FVI(J) = 1.0 
FVIC(J) = 0.0 

RFJ = RFGX*PL(J)*V(J)*V(J)*RH0(ND2) 

RF( J) = DSIGN(RFJ,V( J)) 

Y(NS5+J) =1.0 

RHO(ND2) = P(ND2)*SMW(1)/RGAS/TKS(3) 

Y(NS7+ND2) = RH0(ND2) 

THE INTERFACE IS MOVING WITHIN THE PIPE 

ELSE IF(FVI(J).GE.0.0) THEN 
NFVI = J 

YPRIME (NS 5+ J) = V( J)/PL( J) 

FVIC(J) = 1.0 - FVI(J) 

RFJ = PL( J)*V( J)*V( J)*(FVI( J)*RFGX*RH0(ND1 ) + 

* FVIC ( J ) *RFLX*RHO ( ND2 ) ) 

RF( J) = DSIGN(RFJ,V( J) ) 

THE INTERFACE IS RETURNING UPSTREAM 
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C 

ELSE IF(FVI(J).LT.O.O) THEN 
IF(J.EQ.l) THEN 
FVI(J) = 0.0 
YPRIME ( NS 5+ J ) =0.0 
FVIC(J) = 1.0 

RFJ = PL ( J ) *V ( J ) * V ( J ) *RHOLIQ*RFLX 
RF( J) = DSIGN(RFJ,V(J)) 

ELSE 

FVI(J-I) = 1.0 + FVI(J) 

FVIC(J-l) = 1.0 - FVI(J-l) 

RFJ = PL(J-1)*V(J-1)*V(J-1)*(FVI(J~1)*RFGX* 
* RH0(ND1-1)+FVIC( J-l )*RFLX*RH0(ND1 ) ) 

RF(J-l) = DSIGN(RFJ,V(J-1)) 

Y(NS5+J-1) = FVI(J-I) 

YPRIME (NS5+J-1) = V( J-1)/PL( J-l) 

RHO(NDl) = RHOLIQ 
Y(NS7+ND1) = RHOLIQ 
FVI(J) = 0.0 
FVIC(J) = 1.0 

RFJ = PL ( J ) *V( J ) *V( J ) *RFLX*RHOLIQ 
RF(J)= DSIGN(RFJ,V(J)) 

Y(NS5+J) =0.0 
YPRIME (NS5+J) =0.0 
ENDIF 
ENDIF 


BOTH NODES ARE LIQUID 


ELSE 

YPRIME ( NS5+J ) =0.0 
FVI(J) = 0.0 
FVIC(J) = 1.0 

RFJ = PL(J)*V(J)*V(J) *RFLX*RHOLIQ 
RF( J) = DSIGN(RFJ,V(J)> 

ENDIF 
100 CONTINUE 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


IF THE VELOCITY OF THE FURTHEST DOWNSTREAM PIPE (NP) IS 
NEGATIVE, AN INTERFACE PASSES INTO THE PIPE. THE SYSTEM 
VARIABLE Y(NS10)=FVI2 FOLLOWS THE PROGRESS OF THE INTER- 
FACE UPSTREAM IN THE SAME WAY FVI(J) FOLLOWS THE OXYGEN/ 
HELIUM INTERFACE. NOTE: 

(1) THE INTERFACE THAT ENTERS THE PIPES FROM THE COMBUS- 
TION CHAMBER CAN ONLY GO AS FAR AS THE UPSTREAM NODE 
OF THE FURTHEST DOWNSTREAM PIPE. 

(2) FVI2 IS VIEWED BACKWORD: AS THE INTERFACE GOES UP THE 
PIPE, FVI2 RANGES FROM 0 TO 0.9999. 
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IF( (V(NP) .LE.O. 0) .OR. (RHO(NN) .LE.RGMX) .OR. 

* (FVI2.GT.FVI2MN)) THEN 
IF(FVI2.GT. 1 .0) THEN 

FVI2 = 0.9999 
YPRIME(NSIO) =0.0 
RHO(NN) = RHOCC 
Y(NS7+NN) * RHOCC 

ELSE IF(FVI2.GT.FVI2MN.OR.V(NP) .LT.0.0) THEN 
YPRIME(NSIO) = -V(NP)/PL(NP) 

RHO(NN) = RHOCC 
Y(NS7+NN) = RHOCC 
ELSE IF(FVI(NP).GE.1.0) THEN 
FVI2 = 0.0 
Y(NS10) = 0.0 
YPRIME(NSIO) =0.0 
RHO(NN) = PA*RGAS IN / TKS ( 3 ) / SMW ( 1 ) 

Y(NS7+NN) = RHO(NN) 

ELSE 

FVI2 =0.0 
YPRIME(NSIO) = 0.0 
RHO(NN) = RHOLIQ 
Y(NS7+NN) = RHOLIQ 
ENDIF 
ELSE 

YPRIME(NSIO) =0.0 
ENDIF 

--- FIND PRESSURE AT EACH NODE. 

(1) IF A NODE IS GASEOUS, PRESSURE IS FOUND USING THE 
PERFECT GAS RELATION. 

(2) IF A NODE IS LIQUID, PRESSURE IS FOUND BY CALCU- 
LATING THE PRESSURE LOSSES FROM UPSTREAM PIPES AND 
SUBTRACTING FROM FURTHEST UPSTREAM PRESSURE. 

(3) SOURCE NODE PRESSURES ARE FIXED. PRESSURE AT NODE 
NN (WHICH CONNECTS THE PIPING SYSTEM TO THE COMBUS- 
TION CHAMBER) IS ALSO FIXED. 

P(l) = PS 

XDPIN = 0.5*OMV*OMV*VOL*VOL/(AO2IN*AO2IN*264. 0*264.0* 

* CDIN*CDIN ) 

IF(FVI(NP).LT.1.0) THEN 

IF(FVI2.GT.FVI2MN. OR. V(NP). LT.0.0) THEN 
DPIN = XDPIN / RHOCC 
ELSE 

DPIN = XDPIN / RHOLIQ 
ENDIF 
ELSE 

DPIN = XDPIN / RHOGAS 
ENDIF 


174 


non 


DPIN = DSIGN(DPIN,V(NP>) 

P(NN) = PA + DPIN 
DO 110 J=2,NP 
ND1 = Nl( J) 

IF(RHO(NDl) .LT.RGMX) THEN 

P(ND1) = RHO(NDl) * RGAS * TKS(3) / SMW(l) 

ELSE 

SRF = 0.0 
RHLSM = 0.0 
DO 120 JX-J-l.NP 

SRF = SRE + RF(JX)*AREA(JX) 

RHLSM = RHLSM + RHOLIQ * FVIC(JX) * PL(JX) 

120 CONTINUE 

RHLSM = RHLSM + FVI2*PL(NP)*(RHOCC-RHOLIQ)+ 

* FVI(J)*PL(J)*RH0(ND1) 

DELP1 = P(NDl-l) - P(NN) 

IF(DABS(SRF) .LT. 1 .0D-04) THEN 
P(NDl) = P(NDl-l) 

ELSE 

P(ND1) = P(NDl-l) - RF(J-1)*AREA(J-1)*DELP1/SRF 
ENDIF 

DELP2 = P(ND1) - P(NN) 

DELPP = DELP2/DFL0AT(NN-ND1) 

DO 130 KX = ND1+1,NN-1 
P(KX) = P(KX-l) - DELPP 
130 CONTINUE 

GO TO 140 
ENDIF 
110 CONTINUE 
140 CONTINUE 
C 

C --- THE EQUATIONS FOR DRH(K) AND DV(J) ARE LINEAR FOR THE 
C DERIVATIVES. A MATRIX EQUATION IS SET UP WITH THE COEF- 

C FICIENTS OF DRH AND DV IN MATRIX ALEQ AND THE RHS IN 

C VECTOR BLEQ. IMSL SUBROUTINE LEQT1F IS CALLED TO SOLVE 

C THE MATRIX EQUATION [ALEQ]*{D} = {BLEQ}. VALUES FOR DV 

C AND DRH ARE RETURNED IN BLEQ. 

C 

C CALCULATE THE COEFFICIENTS ALEQ ( NP+K , NP+K ) = SUM1(K) 

C 

C --- RESET SUM1 
C 

DO 150 K=1 ,NN 
150 SUM1(K) = 0.0 

CALCULATE SUM1 

DO 160 J=1,NP 

SUM1(N1(J)) = SUM1(N1(J)> - AREA( J)*V( J) 

SUM1(N2( J)) = SUM1(N2( J)) + AREA(J)*V(J) 
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160 CONTINUE 

DO 161 K-l.NN 

IF(DABS(SUM1(K)) .LT.SM1MIN) SUM1(K)=DSIGN(SM1MIN, 

* SUM1(K) ) 

161 CONTINUE 

— RESET ALEQ AND BLEQ 

NVP = NN + NP 
DO 170 1=1, NVP 
BLEQ(I) = 0.0 
DO 180 J=1 ,NVP 
ALEQ(I, J) = 0.0 
180 CONTINUE 
170 CONTINUE 

--- FILL ALEQ AND BLEQ 

MLIQ = 0 
DO 190 J-l.NP-l 
ND1 = Nl( J) 

ND2 = N2( J) 

BOTH NODES ARE GASEOUS 

IF(RHO(ND2) .LT.RGMX) THEN 

ALEQ(J.J) = 0.5*PL( J)*(RH0(ND1) + RHO(ND2) ) 

ALEQ( J , NP+ND1 ) = 0.5 * PL(J) * V(J) 

ALEQ( J , NP+ND2 ) = 0.5 * PL(J) * V(J) 

BLEQ(J) = P(ND1) - P(ND2) - RF(J) 

ALEQ ( NP+ND 1 , NP+ND 1 ) = SUMl(NDl) 

ALEQ(NP+ND1, J) = - RHO(NDl) * AREA(J) 

ALEQ(NP+ND2, J) = RH0(ND2) * AREA(J) 

--- ONLY THE UPSTREAM NODE IS GASEOUS 

ELSE IF(MLIQ.LT.l) THEN 
MLIQ = 1 
JHE = J 

ALEQ( J, J) = RHLSM 

ALEQ ( J, NP+ND 1) = PL(J) * FVI(J) * V(j) 

ALEQ( J,NP+NN) = PL(NP) * FVI2 * V(NP) 

BLEQ(J) = P(ND1)-AREA(NP)*P(NN)/AREA(J)-SRF/AREA(J)- 

* RHOLIQ*V( J)*V( J)-RHO(NN)*V(NP)*V(NP)*AREA(NP)/ 

* AREA(J) 

ALEQ ( NP+ND 1, NP+ND 1) = SUMl(NDl) 

ALEQ ( NP+ND 1,J) = - RHO(NDl) * AREA(J) 

--- BOTH NODES ARE LIQUID 
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ELSE 

ALEQ(J.J) = 1.0 

ALEQ(J.J-l) = -AREA(J-1)/AREA(J) 

ALEQ ( NP+ND 1 , NP+ND 1 ) =1.0 
ENDIP 
190 CONTINUE 
C 

C — THE FURTHEST DOWNSTREAM PIPE IS HANDLED SEPARATELY BECAUSE 
C OF THE POSSIBILTY OF A BACKFLOW INTERFACE. 

C 

C --- THE UPSTREAM NODE IS LIQUID 
C 

IF(RHO(NN-l) .GT.RGMX) THEN 
ALEQ(NP.NP) =1.0 

ALEQ(NP.NP-l) = -AREA(NP-1)/AREA(NP) 

ALEQ ( NP+NN - 1 , NP+NN - 1 ) = 1.0 
ALEQ ( NP+NN , NP+NN ) = 1.0 

-—THE FURTHEST DOWNSTREAM NODE IS LIQUID 

IF(FVI2.LE.0.0) THEN 
BLEQ( NP+NN) =0.0 

—A BACKFLOW INTERFACE EXISTS 

ELSE 

BLEQ( NP+NN) = RPX 
ENDIF 

BOTH NODES ARE HELIUM 

ELSE IF(FVI(NP).GE.1.0) THEN 

ALEQ(NP.NP) = 0 . 5*PL(NP)*(RHO(NN-l ) + RHO(NN) ) 

ALEQ( NP, NP+NN- 1) = 0.5 * PL(NP) * V(NP) 

ALEQ (NP, NP+NN) = 0.5 * PL(NP) * V(NP) 

BLEQ(NP) = P(NN-l) - P(NN) - RF(NP) 

ALEQ( NP+NN- 1, NP+NN- 1) = SUMl(NN-l) 

ALEQ ( NP+NN -1,NP) = -RHO(NN-l) * AREA(NP) 

ALEQ( NP+NN, NP+NN) =1.0 
BLEQ( NP+NN) = YPRIME(NS3) * RHO(NN) 

AN HELIUM/OXYGEN INTERFACE EXISTS IN THE PIPE 

ELSE IF(FVI2.LE.0.0) THEN 

ALEQ(NP,NP)=PL(NP)*(RHO(NN-l)*FVI(NP) + RHOLIQ*FVIC(NP) ) 
ALEQ( NP, NP+NN- 1) = PL(NP) * V(NP) * FVI(NP) 

BLEQ(NP) = P(NN-l) - P(NN) - RF(NP) 

ALEQ( NP+NN- 1, NP+NN- 1) = SUMl(NN-l) 

ALEQ ( NP+NN -1,NP) = -RHO(NN-l) * AREA(NP) 

ALEQ (NP+NN, NP+NN) = 1.0 
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C --- TWO INTERFACES EXIST IN THE LAST PIPE (BACKFLOW AND HE/OX) 
C 

ELSE 

ALEQ(NP,NP) = PL(NP)*(RHO(NN-l)*FVI(NP)+(FVIC(NP)-FVI2) 

* * RHOLIQ + RHOCC * FVI2) 

ALEQ(NP, NP+NN- 1) = PL(NP) * FVI(NP) * V(NP) 
ALEQ(NP,NP+NN) = PL(NP) * FVI2 * V(NP) 

BLEQ(NP) = P(NN-l) - P(NN) - RF(NP) 

ALEQ ( NP+NN - 1 , NP+NN - 1 ) = SUMl(NN-l) 

ALEQ( NP+NN- 1 ,NP) = -RHO(NN-l) * AREATA) 

ALEQ( NP+NN, NP+NN) ■ 1.0 
BLEQ(NP+NN) = YPRIME(NS3) * RHOCC 
ENDIF 

ALEQ(NP+1,NP+1) =1.0 
ALEQ(NP+1 , 1 ) =0.0 
IERP = 10 

CALL LEQT 1 F ( ALEQ , 1 , NVP , 50 , BLEQ , 0 , WKAREA , IERP ) 
IF(IER.GT.IOO) THEN 

WRITE (30,1300) IERP,TM,TK,PA,RHOCC,FVI2,YPRIME(NS10) 
WRITE (30,1310) ( J,FVI( J) ,YPRIME(NS5+J) ,V( J) ,DV( J) , 

* J=1 ,NP) 

WRITE (30,1320) (K,RHO(K) ,P(K) ,DRH(K) ,K=1 ,NN) 

STOP 

ENDIF 

DO 800 J=1 ,NP 

DV( J) = BLEQ(J) 

YPRIME ( NS 3+ J ) = DV( J) 

800 CONTINUE 

DO 900 K=1,NN 

DRH(K) = BLEQ(NP+K) 

YPRIME ( NS 7+K ) = DRH(K) 

900 CONTINUE 

IF(NDEBUG( 12) ) THEN 

WRITE (30,1400) FVI 2, YPRIME (NS10 ) ,SRF 
WRITE (30,1500) ( J,V( J) ,DV( J) ,RF( J) , J=1 ,NP) 

WRITE (30,1600) (K,RHO(K) ,DRH(K) ,P(K) ,K=1,NN) 

ENDIF 

RETURN 

1000 FORMAT(//2X , ' ***** OUTPUT FROM SUBROUTINE PIPES *****'//) 
1300 FORMAT(//2X, 1 ***** CONVERGENCE OF PIPING SYSTEM FAILED', 

* • ****** //,4X,'LEQT1F RETURNED WITH IER = ',14, 

* ', EXECUTION HALTED 7, 4X, 'TIME = ',1PD11.4, 

* 4X , ' CHAMBER TEMPERATURE = \0PF10.2,’ K'/, 

* 4X, 'CHAMBER PRESSURE = ’,1PD11.4,' PA'/, 

* 4X, 'CHAMBER DENSITY = '.D11.4,' KG/M**3'/, 

* 4X, ' FVI2 = ' ,D11 .4/ , 

* 4X, 'YPRIME (FVI 2) = ’,D11.4) 

1310 FORMAT ( / / , 4X , 4HPI PE , 5X , 3HFVI ,4X,9HD(FVI)/DT,2X, 8HVELOCITY , 

* 5X, 7HD(V)/DT// , (6X, II , 2X, 1PD10. 3,X, 1PD10. 3,X, 1PD10. 3,X, 


* 1PD10.3)) 

1320 FORMAT ( / / 4X , 4HN0DE , 3X , 7HDENSITY , 3X , 8HPRESSURE , 4X , 

* 9HD(RHO)/DT/ / , (6X, II , 2X, 1PD10. 3,X, 1PD10 .3,X,1PD10.3)) 
1400 F0RMAT(/4X' PIPES DERIVATIVE ESTIMATE SUMMARY'//, 

* 4X,' POSITION OF THE BACKFLOW INTERFACE (FVI2) = 

* 1PD11 .4/ ,4X, 'D(FVI2)/DT = ' ,D11.4,4X, 'SRF 

* Dll. 4) 

1500 F0RMAT(/4X, ' PIPE' ,2X, 1 V ',2X,' D(V)/DT ',2X, 

* * RF 1 // , (5X,I2,3X, 1PD11 .4,2X,D11 .4,2X,D11 .4) ) 

1600 FORMAT (/4X, 'NODE' ,2X,' RHO ',2X,' D(RHO)/DT ',2X, 

* ' P '// , (5X,I2,3X, 1PD11 .4,2X,D11 .4,2X,D11 .4)) 

END 
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SUBROUTINE RADSS 
C 

INCLUDE ' INIT .FOR/ LIST' 

DIMENSION Gl(20) ,BLEQ( 10) ,WKAREA( 10) ,ALEQ( 10,10),S2I(10) 

C 

• INCLUDE 'COMMONSS. FOR/LIST ' 

C 

ITER = 0 
C 

SMNS = 0.0 
DO 100 1=1, NS 

• SMNS = SMNS + S2(I) 

100 CONTINUE 

DO 110 I=NS1 , NSRATE 
S2I(I) = S2(I) 

110 CONTINUE 

— - DEFINE THE G1 ARRAY: G1(I) = MASS INFLOW OF SPECIES I FROM 
THE PIPES. 

Gl(2) = FMV / SMW(2) 

Gl(3) = OMV / SMW(3) 

ITERATION BEGINS 

120 CALL RATES 
ICHK = 0 
ITER = ITER + 1 
DO 140 I=NS 1 , NSRATE 

BLEQ(I-NS) = - G1(I) - ARR(I) + EMV * S2(I) 

DO 130 K=NS1, NSRATE 

ALEQ(I-NS.K-NS) = RIP(I,K) 

130 CONTINUE 

ALEQ(I-NS.I-NS) = ALEQ(I-NS,I-NS) - EMV 
140 CONTINUE 

— - LEQT1F SOLVES THE MATRIX EQUATION : [ALEQ] {DS2}={BLEQ} AND 
RETURNS THE ANSWER IN THE VECTOR {BLEQ} . 

CALL LEQT1F(ALEQ, 1 ,NRAD, 10,BLEQ,IDGT,WKAREA, JER) 

INCREMENT ALL S2 BY DS2 AND KEEP THE RESULT WITHIN THE BOUNDS 

( S2MIN < S2(l) < 1.0/SMW(I) ) 

SM = SMNS 

DO 150 I=NS1, NSRATE 
ALG = 1 . 0/SMW( I ) 

SCHK = DABS(BLEQ(I-NS)/S2(I)) 

IF ( SCHK . GT . EPSISS ) ICHK = 1 
S2(I) = S2(l) + RLX*BLEQ(I-NS) 


c-3 
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c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


IF(S2(I).LT.S2MIN) S2(I) = S2MIN 
IF(S2(I) .GT.ALG) S2(I) = ALG 
SM = SM + S2(I) 

150 CONTINUE 

SMINV = 1.0/ SM 

--- STOP THE ITERATION IF ITERATION COUNTER ITER IS GREATER 
THAN ITMAXS. IF THE SPECIES MOLE NUMBERS FALL BELOW THE 
MINIMUM ALLOWABLE MOLE NUMBER THE ITERATION CANNOT CON- 
VERGE. IF THIS TAKES PLACE THE MOLE NUMBERS ARE SET TO 
MINIMUM ALLOWABLE MOLE NUMBER AND EXECUTION IS ALLOWED 
TO CONTINUE. 

IF ( ITER. GT. ITMAXS) THEN 
ISTOP = 0 

S2CHK = S2MIN*100.0 
DO 160 I = NS1 ,NSRATE 

IF(S2(I) .LE.S2CHX) THEN 
ISTOP = 1 
S2(I) = S2MIN 
ENDIF 

160 CONTINUE 

IF (ISTOP. EQ.l) GO TO 170 

WRITE (30,1000) ITMAXS, TM.TK, PA, RHOCC,DLIQ,VR 
WRITE (30,1100) (ANAM(I),S2(I),I=1,NS) 

WRITE (30,1200) (ANAM(I),S2(I),BLEQ(I-NS),I=NS1,NSRATE) 
WRITE (30,1300) 

STOP 

ENDIF 

---IF ICHK - 1 THE CONVERGENCE CRITERION WAS NOT MET 
IF(ICHK) GO TO 120 

--- CALCULATE THE RATE OF PRODUCTION OF THE PRODUCTS AND INFER 

THE RATE OF CONSUMPTION OF THE REACTANTS FROM THE CONSERVATION 
OF ATOMS 


C 

170 RATE ( 4 ) = ARR(A) 

RATE(2) = -RATE(4) 

• RATE ( 3 ) = -RATE(4)*0 . 5 

RATE(l) - 0.0 

C 

C --- OUTPUT FOR NDEBUG(ll) ■ 1 
C 

IF(NDEBUG(11)) THEN 

• WRITE (30,1A00) ITER,RLX,EPSISS 

WRITE (30,1500) TK,PA,RHOCC,OMV,FMV,PMV 
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WRITE (30,1600) 

WRITE (30,1700) (ANAM(l),S2l(l),I=NSl,NSRATE) 

WRITE (30,1800) 

WRITE (30,1900) (ANAM(I) ,S2(I) ,AKR(I) ,RATE(I) ,1=1 ,NSRATE) 
ENDIF 
C 

RETURN 

C 

1000 FORMAT(/2X, ’***** IN RADSS, ITMAX = ',17,' EXCEEDED *****'//, 

* 4X, 'RADICAL MOLE NUMBER ITERATION DID NOT CONVERGE FOR:'//, 

* 4X , ' TIME = MPD11.4,' SEC'/, 

* 4X, 'CC TEMPERATURE = '.0PF10.3,' K'/, 

* 4X, ' CC PRESSURE = '.1PD11.4,' PA'/, 

* 4X, ' CC DENSITY = '.Dll. 4,' KG/M**3'/, 

* 4X, 'CC LIQUID DENSITY = ',D11.4,' KG/M**3’/, 

* 4X, ' CC VAPOR. RATE = ',D11.4,' KG/M**3 SEC'//) 

1100 FORMAT (4X, 'MAJOR SPECIES MOLE NUMBER (KGMOLES/KG) ' // , 

* (8X,A4, 13X,D11 .4) ) 

1200 FORMAT(//4X,' SPECIES MOLE NUMBER (KGMOLES/KG) INCREMENT'//, 

* (5X,A4, 7X,D11 .4, 13X.D11 .4)) 

1300 FORMAT(//2X, ****** EXECUTION HALTED *****') 

1400 FORMAT (/2X, ****** OUTPUT FROM SSRATE *****',//, 

* 4X, 'SSRATE SUCCEEDED AFTER ',14,' ITERATIONS.'/, 

* 4X, 'RELAXATION FACTOR = '.1PD8.2/, 

* 4X, 'PERCENT TOLERANCE = ’,D8.2) 

1500 FORMAT (4X,' TEMPERATURE =',F8.2/, 

* 4X, 'PRESSURE = '.1PD12.4/, 

* 4X, 'DENSITY = '.D12.4/, 

* 4X, 'OMV, FMV, PMV = ’ ,3(D12.4,2X) ) 

1600 F0RMAT(/4X, ’ SPECIES' ,2X, 'INITIAL MOLE NO.’//) 

1700 FORMAT(5X,A4,3X,D12.5) 

1800 FORMAT (/4X, ' CONVERGED SOLUTION: ' // , 

* 4X,' SPECIES ',2X,' EFFECTIVE S2',2X,'SS RXN RATE ',2X, 

* 'EFFECTIVE RXN 

RATE ' / , 1 2X , ' KGMOLE ( I ) / KG ' , X , ’ KGMOLE/M * 3*SEC ' , X , 

* ' KGMOLE/M *3*SEC 7/ ) 

1900 FORMAT(5X,A4,3X,lPD12.5,2X,D12.5,2X,D12.5) 

END 
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SUBROUTINE RATES 


--- CALCULATES RATE OF PRODUCTION OF EACH SPECIES, ARR(I) , 
AND THE DERIVATIVE OF THE RATE OF PRODUCTION OF EACH 
SPECIES, I, WITH RESPECT TO THE MOLE NUMBER OF EACH 
SPECIES, K , RIP(I,K) . 

INCLUDE 'INIT. FOR/ LIST’ 

DIMENSION RFP(20,20),RRP(20,20),RF2(20),RR(20) 

INCLUDE ' COMMONSS . FOR/ LIST ' 

DATA BIG/ 46. 051/ 


--- DETERMINE DENSITY AND OFTEN USED CONSTANTS 

TLN * DLOG(TK) 

TKSQ = TK * TK 

TKSQIN = 1.0/TKSQ 

RHSM = PA * RGASIN * TKINV 

RHOP = RHSM * SMINV 

RHSQ = RHOP * RHOP 

RHINV * 1.0 / RHOP 

--- FIND THE FORWARD AND REVERSE RATES OF EACH OF THE ELEMENTARY 
REACTIONS, J, FOR J=1,JJ (JJ = NO. OF ELEM RXNS) 

DO 10 J=1,JJ 
RF2( J) = 0.0 
RR( J) = 0.0 
I = ID( 1 , J) 

K = ID(2,J) 

M = ID(3, J) 

N = ID(4, J) 

MODE = MODR(J) 

R1 = 0.0 
R2 = 0.0 

TX1 = TACT ( J ) *TKINV - BX(J) - TEN( J)*TLN 
TX2 = TACT2 ( J ) *TKINV - BX2(J) - TEN2( J)*TLN 
R1 = DEXP(-TXl) 

R2 = DEXP( -TX2) 

SECTION FOR FORWARD REACTION RATES 

IF(DABS(TX1) .gt.big) then 

IF(NDEBUG( 13) )WRITE( 30,1000) ITER,J 
GO TO 50 
ENDIF 
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20 


S2(K) 


R1 = R1 * X1(J) 

IF(M0DE-2) 20,30,40 

RF2( J) = R1 * S2( I) * RHSQ * 

GO TO 60 

30 RF2( J) = R1 * RHSM * S2(I) *RHOP 

GO TO 60 

40 RF2( J) = R1 * RHSM * S2(l) * RHOP 

IF(K.NE.O) RF2( J) = RF2(J) * S2(K) * RHOP 
GO TO 60 

50 RF2( J) = 0.0 _ 

60 CONTINUE 

--- SECTION FOR REVERSE RATES 

IF(DABS(TX2) .gt.big) then 

IF(NDEBUG(13)) WRITE(30, 1100) ITER, J 
GO TO 100 
ENDIF 

R2 = X2( J) * R2 
IF(MODE-2) 70,80,90 

70 RR( J) = R2 * S2(M) * RHSQ * S2(N) 

GO TO 10 

80 RR( J) = R2 * RHSM * S2(M) * RHSQ * S2(N) 

GO TO 10 

90 RR( J) = R2 * RHSM * RHOP * S2(M) 

GO TO 10 

100 RR(J) = 0.0 

10 CONTINUE 

--- FIND THE RATE OF PRODUCTION OF EACH SPECIES I, ARR(I) . 

DO 110 1=1 .NSRATE 
ARR(I) = 0.0 
DO 120 J=1,JJ 

ARR(I) = ARR(I)+DFL0AT(ISIDE(I,J))*(RF2(J)-RR(J)) 

120 CONTINUE 
110 CONTINUE 

--- FIND THE DERIVATIVES OF THE FWD AND REV RXN RATES OF EACH 
ELEMENTARY RXN WITH RESPECT TO EACH MOLE NO. (I.E. 

RFP( J,K)=D(RF2( J) )/D(S2(K) ) AND RRP(J,K)=D(RR(J))/D(S2(K))) 

DO 130 J=1,JJ 
DO 140 K=1,NSRATE 
RFP( J,K) = 0.0 
RRP( J,K) = 0.0 
IF(ISIDE(K,J).LT.O) 
RFP(J,K)=-DFLOAT(ISIDE(K,J))*RF2(J)/S2(K) 

IF( ISIDE(K, J) .GT.O) RRP(J,K)= DFLOAT(ISIDE(K, J) )* 
RR(J)/S2(K) 


n n 


140 CONTINUE 
130 CONTINUE 

— - FIND THE DERIVATIVE OF THE RATE OF PRODUCTION OF EACH SPECIES 
WITH 

C RESPECT TO THE MOLE NO. OF EACH SPECIES 
(RIP(I,K)=D(ARR(I))/D(S2(K) ) 

c 

DO 150 1=1 ,NSRATE 
DO 160 K=1,NSRATE 
RIP(I,K) = 0.0 
DO 170 J=1,JJ 

RIP(I,K) = RIP(l,K)+DFLOAT(lSIDE(I,J))*(RFP(J,K)-RRP(J,K)) 
170 CONTINUE 
160 CONTINUE 
150 CONTINUE 
RETURN 

1000 FORMAT (4X, ' IN RATES, FOR ITERATION ' ,14,'FWD RATE = 0 FOR J = 

t 

* 12) 

1100 FORMAT (4X, 1 IN RATES, FOR ITERATION ’,14,’REV RATE = 0 FOR J = 


nno nnoo non nnnno o; oonnnoo 


SUBROUTINE YCAL( Y, WK, IWK,N,NWK) 


TSTR2 RUNS THE CALLING OF THE INTEGRATION ROUTINE AND 
THE OUTPUTTING OF DATA TO THE OUTPUT AND PLOT FILES 

INCLUDE 'INIT. FOR/LIST' 

REAL*4 SDUMMY(4) 

EXTERNAL CSLOPE 

— . DIMENSION Y(N) , IWK(N) ,WK(NWK) , WK2(4, 100) ,S2P( 10) ,CHG(50) 

INCLUDE ' COMMONSS . FOR/LIST ' 

COMMON / GEAR / DUMMY(48),SDUMMY,IDUMMY(38) 

---PUT MOLE NUMBERS, TEMPERATURE AND DENSITY INTO LOG VARIABLE 
FORM AND INSERT THEM INTO THE VECTOR Y. DENSITY OF LIQUID IN 
THE COMBUSTION CHAMBER IS LEFT IN NON-LOG FORM. 

DO 100 1=1, NS 
100 Y(I) = DLOG(S2(I)) 

Y(NS1) - XMASS/VOL 
Y(NS2) = DLOG(TK) 

Y(NS3) - DLOG(PA) 

VH2 = EMS ( 2 ) / AH2IN/ RHOH2/ 264.0 

PLACE THE PIPING SYSTEM VARIABLES INTO THE VECTOR Y. 

DO 110 1=1, NP 

Y(NS3+I) = V(I) 

Y(NS5+I) =0.0 

IF(N2(I) .EQ.NN) OMV = OMV + RHO(NN)*AREA(I)*V(I) 

110 CONTINUE 

DO 120 1=1, NN 

Y(NS7+I) = RHO(I) 

120 CONTINUE 

--- GET GOOD INTITIAL ESTIMATES FOR ATOM AND RADICAL MOLE 
NUMBERS VIA A CALL TO RADSS 

IF(KINET.EQ.l) THEN 
CALL RADSS 
ENDIF 

BEGIN THE INTEGRATION LOOP 

NSTEPS = JIDINT((TMEND-TMI)/TMPRNT) 

DO 150 NX = 1, NSTEPS 

TMNXT = TMNXT + TMPRNT 
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IF ( TMNXT . GT . TMEND ) TMNXT = TMEND 
IF(NDEBUG(5)) WRITE (30,1000) NX, NSTEPS,TM, TMNXT 
IF(NDGEAR) THEN 

CALL DTSTR( N , CSLOPE , TM , STEP , Y , TMNXT , WK2 , NDEBUG ) 

ELSE 

CALL DGEAR (N, CSLOPE, CHGJ.TM, STEP, Y, TMNXT, EPSI.METH, 

* MITER, INDEX, IWK,WK,IER) 

IF(IER.GE. 128) THEN 

WRITE (30,1100) IER.TM 
DO 140 1=1, NV 

140 — CHG(I) = WK(IDUMMY(11)+I)/WK(I) 

WRITE (30,1200) (I,Y(I),CHG(I),I=1,NV) 

CALL OUTPT(Y,TM,N,NSTEPS) 

STOP 
END IF 
END IF 

IF (NDEBUG(14)) 

* WRITE (6,1300) TM,STEP,TK,PA,RH0CC,(V(J),J=1,3), 

* (RHO(K) ,K=1 ,3) , (FVI( J) , J=NP-2,NP) ,V(NP) ,FVI2, 

* (S2(I) ,1=1 ,4) 

ON RETURN FROM THE INTEGRATION ROUTINE OUTPUT DEPENDENT 

VARIABLES AND SIGNIFICANT PROBLEM PARAMETERS TO THE TEXT 
AND PLOT OUTPUT FILES. 

CALL OUTPT(Y,TM,N,NSTEPS) 

IF(TM.GT.TMPER.AND.IPER.LT.l) THEN 
IPER = 1 
DMCX = DMC 
DMC = XPER * DMC 

ELSE IF ( TM . GT . TMPER . AND . IPER . LE . 5 ) THEN 
IPER = IPER + 1 
ELSE IF(IPER.GT.5) THEN 
DMC = DMCX 
ENDIF 

150 CONTINUE 
C 

RETURN 

1000 FORMAT(/2X, '***** INTEGRATING ROUTINE CALLED *****'//, 

* 4X , ' CALL FOR STEP ’,17,’ OUT OF ',112,' STEPS’/, 

* 4X , ' STEP BEGINS AT TIME = ',1PD12.5, 

* ' AND ENDS AT TIME = '.1PD12.5) 

1100 FORMAT (2X, '***** ERROR CONDITION ON RETURN FROM DGEAR 
*****' / / 

* 4X, 'EXECUTION STOPPED. IER = ',I4,X,'TIME = ',1PD12.5//, 

* 2X,' VARIABLE’, 2X,' VALUE AT HALT' ,2X, 'RELATIVE CHANGE'//) 

1200 FORMAT(5X,I2,5X, 1PD12.5,4X,D12.5) 

1300 FORMAT(/2X, 'ON RETURN TO YCAL: TIME = ',1PD8.2,' STEP = 

' ,D8.2/, 

* 2X , ' TEMP = ' ,0PF10.2,X, 'PRESS = ' , 1PD8.2,X, 'DENS = \D8.2/, 


t 


2X, ' V( 1 : 3) = ' , 3D9.2/, 
2X, 1 RHO( 1:3) = \3D9.2/ } 
2X,'FVI(NP-2:NP) = ' ,3D9.2/, 
2X,'V(NP) = ' ,D9.2,2X, 'FVI2 
2X, ' S2 ( 1 :4) = 1 ,4D9.2) 

END 




188 




' ,D9.2/ , 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 


C 


SUBROUTINE THERPS (NS, XTK) 

THERPS CALCULATES THE SPECIFIC MOLAR INTERNAL ENERGY, 
ENTHALPY AND SPECIFIC HEAT AT CONSTANT PRESSURE FOR EACH 
SPECIES, 1=1, NS. 

H(I) = ENTHALPY ( J/KGMOLE) 

U(I) = INTERNAL ENERGY (J/KGMOLE) 

CP(I) = SPECIFIC HEAT AT CONSTANT PRESSURE (J/KGMOLE K) 

THERPS IS MODELLED AFTER SUBROUTINE HCPS USED BY DR. PAUL 
GEORGE II (PHD DISERTATION , PURDUE UNIVERSITY, 1982). 

IMPLICIT REAL* 8 (A-H.P-Z) 

IMPLICIT INTEGER** (I-O) 

DIMENSION H1(20),H2(20),H01(20),H02(20),CP1(20),CP2(20), 

* CP01(20),CP02(20),S1(20),S2(20),S01(20),S02(20),U1(20), 

* U2(20) 

COMMON 

* / CPROP / Z(200) ,H( 10) ,CP( 10),U(10),S(10) ,DCP( 10), HO (10), 

* CP0(10),S0(10) 


DATA ICIT/14/.RGAS/8314.4/ 

IT=0 

IF (XTK.LT. 950.0) THEN 
IT-7 

ELSE IF(XTK.GT. 1050.0) THEN 
IT = 0 
ELSE 

GO TO 70 
ENDIF 

XTKSQ=XTK**2 
XTKCU=XTKSQ*XTK 
XTK4=XTKCU*XTK 
XTKINV=1 . /XTK 
XTLN=DLOG ( XTK ) 

RTK=RGAS*XTK 

H(I) , U(I), CP(I), DCP(I) REQUIRED 

DO 60 1=1, NS 

K=IT+ICIT*( I- 1 ) 

CP1X=Z(K+1) 

CP2X=XTK*Z(K+2) 

• CP3X=XTKSQ*Z ( K+3 ) 

CP4X=XTKCU*Z ( K+4 ) 
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CP5X=XTK4*Z(K+5) 

HO(I)=Q . 2*CP5X+0 . 25*CP4X+Q . 33333333*CP3X+0 . 5*CP2X+CP1X 

* +XTKINV*Z(K+6) 

S0(I)=0T25*CP5X+0 . 33333333*CP4X+0 . 5*CP3X+CP2X+CP1X*XTLN 
CPO(I) = CP1X + CP2X + CP3X + CP4X + CP5X 
H(I) = HO(I) * RTK 
U(I) = H(I) - RTK 
S(I) = SO(I) * RGAS 
CP(I) = CP0(I)*RGAS 
60 CONTINUE 

RETURN 

70 XTKSQ=XTK**2 
XTKCU=XTKSQ*XTK 
XTK4=XTKCU*XTK 
XTKINV=1 . /XTK 
XTLN=DLOG ( XTK ) 

RTK=RGAS*XTK 

H(I), U(I), CPU), DCP(I) REQUIRED 

DO 80 ITK = 1,2 
IT - 0 

IF(ITK.EQ.2) IT * 7 
DO 90 1=1, NS 

K=IT+ICIT* ( I - 1 ) 

CP1X=Z(K+1) 

CP2X=XTK*Z ( K+2 ) 

CP3X=XTKSQ*Z(K+3) 

CP4X=XTKCU*Z ( K+4 ) 

CP5X=XTK4*Z ( K+5 ) 

H0( I )=0 . 2*CP5X+0 . 25*CP4X+0 . 33333333*CP3X+0 . 5*CP2X+CP1X 

* +XTKINV*Z ( K+6 ) 

S0( I)=0 . 25*CP5X+0 . 33333333*CP4X+0 . 5*CP3X+CP2X+CP1X*XTLN 
CPO(I) = CP1X + CP2X + CP3X + CP4X + CP5X 
H(I) = H0( I ) * RTK 
U(I) = H(I) - RTK 
S(I) = S0(I) * RGAS 
CP(I) = CP0(I)*RGAS 
IF(ITK.EQ.l) THEN 
HO 1 ( I ) = H0(I) 

S01(I) = S0(I) 

CPOl(I) = CPO(I) 

H1(I) = H( I) 

CP1(I) = CP(I) 

U1(I) = UU) 

S1(I) = S(I) 

ELSE 

H02(I) = HO(I) 



S02 ( I ) = SO(I) 

CP02(I) = CPO(I) 

H2( I) = HU) 

CP2(I) = CP(I) 

U2(I) = U(I) 

S2(I) = S(I) 

ENDIF 
90 CONTINUE 
80 CONTINUE 

FR1 = (XTK - 950.0)/100.Q 
FR2 = (1050.0 - XTK)/100. 0 
DO 100 1=1, NS 

H0(I) = FR1*H01(I) + FR2*H02(I) 
S0(I) = FR1*S01(I) + FR2*S02(I) 
CPO(I) = FR1*CP01(I) + FR2*CP02(I) 
H(I) = FR1*H1(I) + FR2*H2(I) 

U(I) = FR1*U1(I) + FR2*U2(I) 

CP(I) = FR1*CP1(I) + FR2*CP2(I) 

100 CONTINUE 


RETURN 


» 


INPUT FILE TSTR.IN 

&RUNNO NMON=06,NDAY=11,NYR=1987,NRUN=1 &END 

&PARAM V0L=0 . 013 , Q= -5 . 5D+07 , S2MIN= 1 . OD- 18 , EMS=0 . 0 , 

18. 00, 0.0,TKS=120. 0,120.0, 200. 0,TMSD=0.0OE+O3, 

RGMX=65 . 5 ,RH0GAS=12 . 4435 ,RH0LIQ=990 . 70 ,RH0H2=9 . 86784, 
AH2IN=1 . 329D-05 ,A02IN=1 . 09682D-05 ,FVI2MN=1 . 0D-08 , 
DLIQMN= 1 . 0D-07 ,DMMIN= 1 . 5D-06 ,DMC=0 . 585 , CFTP=0 . 96738 , 
PHG=1. 03393D+06 ,TMPER=5 .0D+02 ,XPER=20. 0,CVR=1 . 23045D-02 
&END 

&INITC XMASS= 1 . 000- 15 ,TK= 5 50 . 00 , PA= 5 . 16964D+06 , 

DM=1.0D-07 &END 

&INDX1 NDEBUG=0, 0,0, 0,0,0, 1,0, 0,0, 0,0, 0,0, 0,0, 0,0,0,0, 
MCON=1,KINET=0 &END 

&ICONS MEIH=2 ,MITER=2 ,NDGEAR=0 , IER= 10 , STEP= 1 . OD- 12 , 

TMI=0 . 0 ,TMPRNT=5 . OD-04 , TMEND= 1 . 0D-01 , INDEX= 1 , 
EPSI=1.0D-05 &END 

&SSRRT EPSISS=1.0D-06,I1MAXS=25,IDGT=0,JER=10 &END 

&INPIP V=9*0. 001, FVI=5*1.0, 4*0. 0,RH0=6*12. 4435, 

4*990 . 70 ,FVI2=0 . 0 , PL= 5*0 . 30 , 0 . 005 , 3*0 . 040 , 

AREA=5*4 . 29085D-05 ,4*2 .02683D-03 ,RFLIQ=0 . 130294, 

RFGAS=0 . 1827 25 ,RFLH0=0 . 836663 , RFGHO= 1 . 14483 , 
RFLV=1.72000D+04,RFGV=1.82725D+04,PS=5.16964D+06, 

CDIN=0 . 50 , SM1MIN=0 . 0D-08 &END 

&INDX2 Nl=l,2,3,4,5,6,7 ,8,9,N2=2,3,4,5,6,7 ,8,9,10, 
NP=9,NN=10,NHEO=5,NFVI=6 &END 


REAC 

HE 

4.0026 

3 

l.OOOOOD+OO 

O.OOOOOD+OO 

G 

O.OOOOOD+OO -3.5416D+06 

REAC 

H2 

2.0160 

2 

l.OOOOOD+OO 9 . 99999D-01 

G 

O.OOOOOD+OO -3.1432D+06 

REAC 

02 

32.0000 

1 

l.OOOOOD+OO 

1.00000D-09 

L 

2 . 00000D+O5 -2.7300D+05 

PROD 

H20 

18.0160 

0 

O.OOOOOD+OO 

1.00000D-09 

G 

O.OOOOOD+OO 

O.OOOOD+OO 

RADI 

H 

1.0080 

0 

O.OOOOOD+OO 

O.OOOOOD-OO 

G 

O.OOOOOD+OO 

O.OOOOD+OO 

RADI 

0 

16.0000 

0 

O.OOOOOD+OO 

O.OOOOOD-08 

G 

O.OOOOOD+OO O.OOOOD+OO 

RADI 

OH 

17.0080 

0 

O.OOOOOD+OO 

0.00000D-08 

G 

O.OOOOOD+OO 

O.OOOOD+OO 

HE 

0 . 25000000D+01 0, 

L 5/66HE 1.00 
.OOOOOOOOD+OO ( 

0.00 0.00 O.G 
3. OOOOOOOOD+OO 

300.000 5000.000 

O.OOOOOOOOD+OO O.OOOOOOOOD+OO 


-0.74537498D+03 0. 915348880+00 0.25000000D+01 0 . OOOOOOOOD+OO O.OOOOOOOOD+OO 3 

O.OOOOOOOOD+OO O.OOOOOOOOD+OO-O. 7453 74980+03 0.915348840+00 4 

0 2.7080 37.0000 298.0 1200.0 6000.0 HE 5 

AR L 5/66AR 1 G 300.000 5000.000 

0 . 25000000E+01 0.0 0.0 0.0 0.0 

-0.74537502E+03 0.43660006E+01 0. 25000000E+01 0.0 0.0 

0.0 0.0 -0. 74537498E+03 0.43660006E+01 

0 3.4180 124.000 300.0 1200.0 5000.0 

CH4 J 3/61C 1H 400 000 OG 300.000 5000.000 1 

0.15027072E 01 0.10416798E-01-0.39181522E-05 0.67777899E-09-0.44283706E-13 2 

-0.99787078E 04 0.10707143E 02 0.38261932E 01-0.39794581E-02 0.24558340E-04 3 

-0.22732926E-07 0.69626957E-11-0.10144950E 05 0.86690073E 00 4 

0. 3.8220 137.0000 

CH3 J 6/69C 1H 3 0 OG 300.000 5000.000 1 
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2.84003270E+00 6.08690860E-03-2 .17403380E-06 3.60425760E-10-2.27253000E-14 2 

1. 64498 130E+04 5.50567510E+00 3.46663500E+00 3.83018450E-03 1.01168020E-06 3 

-1.88592360E-09 6.68031820E-13 1.63131040E+04 2 .41721920E+00 4 

CHO C 1H 10 1 OG 300.000 5000.000 1 

3 . 61363839E+00 3.26359808E-03-1.29355163E-06 2.38238898E-10-1.64668247E-14 2 

3. 1811314 1E+03 5.2216507 8E+00 4.07467474E+00-2.20831318E-03 1. 12814616E-05 3 

- 1.07 609 505E-08 3.39297805E-12 3.33218275E+03 4.03515416E+00 4 

CH20 J 3/61C 1H 20 1 OG 300.000 5000.000 1 

2 . 83642490E+00 6.86052980E-03-2.68826470E-06 4.79712580E-10-3.21184060E-14 2 

-1.52360310E+04 7.85311690E+00 3.79637830E+00-2.57017850E-03 1.85488150E-05 3 

-1.78691770E-08 5.55044510E-12-1.50889470E+04 4.75481630E+00 4 

CH30 C 1H 30 1 OG 300.000 5000.000 1 

5.23759870E+00 7.99608790E-03-3.96765581E-06 8. 31862800E-10-6.26533814E-14 2 

1 . 09662525E+03- 5 . 89948723E+00 2 .04832571E+O0 1.15498951E-02-3.75331995E-06 3 

2 . 54155659E-13 1.69432993E-24 2.63293614E+03 1.27359509E+01 4 


CO J 9/65C 10 100 000 OG 300.000 5000.000 

0.29840696E 01 0.14891390E-02-0.57899684E-06 0. 10364577E-09-0.69353550E-14 
-0.14245228E 05 0.63479156E 01 0.37100928E 01-0. 16190964E-02 0.36923594E-05 
-0.20319674E-08 0.23953344E-12-0.14356310E 05 0.29555351E 01 

0 3.5900 110.0000 298.0 1200.0 6000.0 CO 

C02 J 9/65C 10 200 000 OG 300.000 5000.000 

0.4460804 IE 01 0.30981719E-02-0.12392571E-05 0.22741325E-09-0. 15525954E-13 
-0.48961442E 05-0.98635982E 00 0.24007797E 01 0.87350957E-02-0. 66070878E-05 
0.2002 186 IE-08 0.63274039E-15-0.48377527E 05 0.96951457E 01 

0 3.9960 190.0000 298.0 1200.0 6000.0 C02 

C02- J12/66C 1.0 2.E 1.00 O.G 300.000 5000.000 

0 . 45454636E+01 0 . 26054317E-02 -0 . 109287 34E-05 0 . 20454421E-09-0 . 14184 542E- 13 
-0.54761969E+05 0. 18317366E+01 0. 347437 38E+01 0. 16913805E-02 0.73533802E-05 
-0.99554249E-08 0.36846715E-11-0.54249051E+05 0.83834333E+01 
0 . 0 . 0 . 

N J 3/77N 1. 0. 0. O.G 300.000 5000.000 

0.24370813E+01 0.13233886E-03-0.90907747E-07 0.22864058E-10-0.13762291E-14 
0 . 56128586E+05 0.45211115E+01 0.25000000E+01-0.31078153E-08 0.83216099E-11 
-0.94278291E-14 0.38108037E-17 0. 56106977E+05 0.41806431E+01 
0. 3. 71. 

N02 J 9/64N 1.0 2.00 0.00 O.G 300.000 5000.000 

0.46240768E+01 0.25260332E-02-0. 10609501E-05 0.19879239E-09-0.13799383E-13 
0 . 22899900E+04 0.13324137E+01 0.34589233E+O1 0.20647063E-02 0.66866069E-05 
-0. 955567 36E -08 0.36195881E-11 0. 28152266E+04 0.83116980E+01 
0 . 0 . 0 . 

HCN L12/69H l.C l.N 1.0 O.G 300.000 5000.000 

0.37068119E+01 0.33382804E-02-0.11913316E-05 0.19992917E-09-0.12826451E-13 
0. 14962637E+05 0.20794907E+01 0.24513559E-K)1 0.87208375E-02-0. 10094203E-04 
0. 672557 12E-08 -0. 17626959E-11 0. 15213004E+05 0.80830088E+01 
0 . 0 . 0 . 


1 

2 

3 

4 
1 
1 
2 

3 

4 
1 
1 
2 

3 

4 

5 
1 
2 

3 

4 

5 
1 
2 

3 

4 

5 
1 
2 

3 

4 

5 
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NH3 J 9/65N l.H 3.00 0.00 0.G 300.000 5000.000 1 

0.24165173E+01 0.61871223E-02-0.21785136E-05 0.37599079E-09-0.24448857E-13 2 

-0.64747187E+04 0. 77043486E+01 0.35912771E+01 0.49388665E-03 0.83449322E-05 ... 3 

-0.83833385E-08 0.27299092E-11-0.66717148E+04 0.22520962E+01 4 

0. 0. 0. 5 

C2H6 C 2H 6 0 OG 300.000 5000.000 1 

1.49270752E+00 1.90694888E-02-7.11984782E-06 1.19069803E-09-7.37266162E-14 2 

-1.13918457E+04 1.36588993E+01 3.06380013E+00 8.59581465E-03 1.64679051E-05 3 

-2.00919579E-08 6.52385805E-12-1. 15878959E+04 6.93085386E+O0 4 

0 4.4180 230.0000 5 


C2H4 C 2H 4 G 300.000 5000.000 

0.34552152E 01 0.11491803E-01-0.43651750E-05 0. 76155095E-09-0.50123200E-13 
0.44773119E 04 0.26987959E 01 0.14256821E 01 0.11383140E-01 0.79890006E-05 
-0.16253679E-07 0.67491256E-11 0.53370755E 04 0.14621819E 02 
0 4.2320 205.0000 

FUEL OIL (LIQUID) C 1H 1.8 G 300.000 5000.000 

0.27797000E 01 

-0.37392400E 04-0.15651223E 02 0.27797000E 01 

-0.37392400E 04-0. 15651223E 02 

NO TRANSPORT 

C3H8 C 3H 8 0 OG 300.000 5000.000 1 

3.07371055E+G0 2.52859519E-02-9.12630575E-06 1.46940810E-09-8.73487271E-14 2 

-1.45316354E+04 7 .48445307E+00 3.21459483E+00 1.63977335E-02 1.68421811E-05 3 

-2. 161117 37E-08 5.77092196E-12-1.42864159E+04 8.64456968E+00 4 

0 5.0610 254.0000 5 

C6H6 C 6H 6 0 OG 300.000 5000.000 1 

9 . 66494563E+00 2.12364775E-02-7.09034255E-06 1.01513328E-09-5.14055768E-14 2 

5 . 17497547E+03-3 .07774861E+01 3.01886004E+00 1.22221561E-02 6.87572729E-05 3 

-9. 9338647 4E-08 4.01152941E-11 8. 100803 15E+03 9.63207440E+00 4 

0 5.27 440. 

C8H18 C 8H 18 0 OG 300.000 5000.000 1 

1 . 14954126E+01 6.35313436E-02-2.71361101E-05 4.80764212E-09-3.02563294E-13 2 

-3 . 35887669E+04-3 . 304295 13E+01 2.27071846E+01-1.27549696E-02 1.04233714E-04 3 

-6 .08467720E-08-9 . 43552892E- 13-3 . 39058757E+04-7 . 78443415E+01 4 

7.451 320. 

H J 3/77H 1 0 0 OG 300.000 5000.000 1 

0 . 25000000E+01 O.OOOOOOOOE 00 O.OOOOOOOOE 00 O.OOOOOOOOE 00 O.OOOOOOOOE 00 2 

0 . 25474390E+05-0 . 4598984 1E+00 0. 25000000E+01 O.OOOOOOOOE 00 O.OOOOOOOOE 00 3 

O.OOOOOOOOE 00 O.OOOOOOOOE 00 0 . 2 5474390E+05-0 . 4598984 1E+00 4 

0 2.7080 37.0000 298.0 1200.0 6000.0 H 1 

H2 J 3/61H 2 0 0 OG 300.000 5000.000 1 

0.31001883E 01 0. 51119458E-03 0.52644204E-07-0.34909964E-10 0. 36945341E-14 2 

-0.87738013E 03-0. 19629412E 01 0.30574446E 01 0.26765198E-02-0.58099149E-05 3 

0.55210343E-08-0. 18122726E-11-0.98890430E 03-0.22997066E 01 4 

0 2.9150 38.0000 300.0 1200.0 5000.0 H2 1 

H20 J 3/61H 20 100 000 OG 300.000 5000.000 1 

0.27167633E 01 0.29451374E-02-0.80224374E-06 0.10226682E-09-0.48472145E-14 2 

-0.29905826E 05 0.66305671E 01 0.40701275E 01-0. 11084499E-02 0.41521180E-05 3 

-0.29637404E-08 0.80702103E-12-0.30279722E 05-0.32270046E 00 4 
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1.2000 2.7100 506.0000 298.0 1200.0 6000.0 H20 

H02 J 3/64H 10 200 000 0G 300.000 5000.000 

0.37866280E 01 0.27885404E-02-0.10168708E-05 0. 17183946E-09-0. 11021852E-13 
0. 11888500E 04 0.48147611E 01 0.35094850E 01 0.11499670E-02 0. 58784259E-05 
-0.77795519E-08 0.29607883E-11 0.13803331E 04 0.68276325E 01 

0 3.4700 107.000 298.0 1200.0 6000.0 H02E 

H202 L 2/69H 20 20 00 OG 300.000 5000.000 

0.45731667E 01 0.43361363E-02-0.14746888E-05 0.23489037E-09-0. 14316536E-13 
-0.18006961E 05 0.50113696E 00 0.33887536E 01 0.65692260E-02-0.14850126E-06 
-0 . 46258055E-08 0.24715147E-11-0.17663147E 05 0.67853631E 01 

UNKNOWN 

N J 3/77N 1 0 0 OG 300.000 5000.000 

0 . 24370811E+01 0. 13233886E-03-0 . 90907754E-07 0 . 22864054E-10-0 . 1376229 IE- 14 
0.56128585E+O5 0.45211111E+01 0.25000004E+01-0.31078154E-08 0.83216097E-11 
-0.94278278E-14 0.38108039E-17 0.56106975E+05 0.41806431E+01 

0. 3.2980 71.4000 300.0 1200.0 5000.0 N 

NO J 6/63N 10 100 000 OG 300.000 5000.000 

0.31890000E 01 0.13382281E-02-0.52899318E-06 0.95919332E-10-0.64847932E-14 
0.98283290E 04 0.67458126E 01 0.40459521E 01-0.34181783E-02 0. 79819 190E-05 
-0. 611393 16E-08 0. 15919076E-11 0.97453934E 04 0.29974988E 01 

0 3.4700 119.0000 298.0 1200.0 6000.0 NO 

N2 J 3/77N 2 0 0 OG 300.000 5000.000 

0 . 28532899E+01 0.16022128E-02-0.62936893E-06 0.11441022E-09-0.78057465E-14 
-0. 8900809 3E+03 0.63964897E+01 0.37044177E+01-0. 14218753E-02 0.28670392E-05 
-0 . 12028885E-08-0 . 13954677E- 13-0 . 10640795E+04 0.22336285E+O1 

0 3.6810 91.5000 300.0 1200.0 5000.0 N2 

0 J 3/770 1 0 0 OG 300.000 5000.000 

0 . Z5342961E+01-0 . 12478170E-04-0 . 125627 24E-07 0.69029862E-11-0.63797095E-15 
0.29231108E+05 0. 4962859 1E+01 0.30309401E+01-0.2Z5Z5853E-02 0.39824540E-05 
-0 . 32604921E-08 0. 10152035E-11 0.29136526E+05 0.26099342E+01 

0 3.0500 106.7000 298.0 1200.0 6000.0 O 

OH J12/700 1H 10 00 OG 300.000 5000.000 

0.29131230E+01 0.95418248E-03-0. 19084325E-06 0. 12730795E-10 0.24803941E-15 
0 . 39647 060E+04 0.54288735E+01 0.38365518E+01-0.10702014E-02 0. 948497 57E-06 
0 . 20843575E-09-0 . Z3384Z65E- 12 0. 36715807E+04 0.49805456E+00 

.7540 2.6620 506.0000 300.0 1200.0 5000.0 OH 

02 J 9/650 2 0 0 OG 300.000 5000.000 

0.36219521E 01 0.73618256E-03-0.19652219E-06 0. 36201556E-10-0.28945623E-14 
-0.12019822E 04 0.36150942E 01 0.36255980E 01-0.18782183E-02 0.70554543E-05 
-0 . 67 63507 IE-08 0.21555977E-11-0.10475225E 04 0.43052769E 01 

0 3.4670 106.7000 300.0 1200.0 5000.0 02 

H20(L) J11/65H 20 100 000 OL 273.150 1000.000 

0.0 0.0 0.0 0.0 0.0 

0.0 0.0 0. 12712782E 02-0.17662790E-01-0.22556661E-04 

0.20820908E-06-0.24078614E-09-0.37483200E 05-0.59115345E 02 

NO TRANSPORT 


1 

1 

2 

3 

4 

1 

1 

2 

3 

4 

1 

2 

3 

4 
1 
1 
2 

3 

4 
1 
1 
2 

3 

4 
1 
1 
2 

3 

4 
1 
1 
2 

3 

4 
1 
1 
2 

3 

4 
1 
1 
2 

3 

4 

5 


C2H6 


02 


C2H4 HZ 


14.41 0.0 Z5000. GLOBAL 
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C3H8 

02 


C2H4 

H2 


14.41 

0.0 

25000. GLOBAL 

C8H18 

02 


C2H4 

H2 


14.41 

0.0 

25000. GLOBAL 

C2H4 

02 


CO 

H2 


12.57 

0.0 

25200. GLOBAL 

C6H6 

02 


CO 

H2 


7.720 

1.0 

19650. GLOBAL 

FUEL OIL02 


CO 

H2 


7.720 

1.0 

19650. GLOBAL 

H20(L) 



H20 



13.90 

0.0 

25900. SUUBERG 

02 

H2 


OH 

OH 


11.903 

0.0 

22661. PRATI 

OH 

H2 


H20 

H 


11.439 

0.0 

5187. PRATT 

02 

H 


OH 

0 


11.677 

0.0 

8712. PRATT 

0 

H2 


OH 

H 


9.352 

0.0 

3903. PRATT 

0 

H20 


OH 

OH 


11.095 

0.0 

9115. PRATT 

H 

H 

M 

H2 


M 

12.699 

-1.15 

0.0 PRATT 

0 

0 

M 

02 


M 

9.672 

-0.278 

0000 PRATT 

0 

H 

M 

OH 


M 

10.627 

-0.0 

-1400. PRATT 

H 

OH 

M 

H20 


M 

10.778 

-0.00 

-252.0 PRATT 

H 

02 

M 

H02 


M 

9.186 

0.0 

503.5 WALD7 

CO 

OH 


H 

C02 


8.336 

0.0 

300. PRATT 

CO 

0 

M 

C02 


M 

8.000 

0.0 

1259. PRATT 

C02 

0 


CO 

02 


10.279 

0.0 

27268. PRATT 

CH4 

0 


CH3 

OH 


13.300 

0.0 

9.000 CGS 

CH4 

H 


CH3 

H2 


14.1 

0.0 

11.900 CGS 

CH4 

OH 


CH3 

H20 


13.500 

0.0 

5.000 CGS 

CH3 

CHO 


CH4 

CO 


11.5 

0.50 

0.0 CGS 

CH3 

CH20 


CHO 

CH4 


10.00 

0.5 

6.000 CGS 

CH3 

0 


CH20 

H 


14.41 

0.0 

2.000 CGS 

CH3 

OH 


CH20 

H2 


12.6 

0.0 

0.0 CGS 

CH3 

H 

M 

CH4 


N 

26.86 

-3.0 

0.0 CGS 

CH3 

02 


CH30 

0 


13.380 

0.0 

28.810 CGS 

CH30 


M 

CH20 

H 

M 

13.700 

0.0 

21.000 CGS 

CH20 


M 

CHO 

H 

M 

16.700 

0.0 

72.000 CGS 

CH20 

0 


CHO 

OH 


13.700 

0.0 

4.600 CGS 

CH20 

H 


CHO 

H2 


13.130 

0.0 

3.760 CGS 

CH20 

OH 


CHO 

H20 


10.5 

1.0 

0.0 CGS 

CHO 


M 

CO 

H 

M 

14.160 

0.0 

19.000 CGS 

CHO 

H 


CO 

H2 


14.300 

0.0 

0.0 CGS 

CHO 

OH 


CO 

H20 


14.000 

0.0 

0.0 CGS 

CHO 

0 


CO 

OH 


14.000 

0.0 

0.0 CGS 


m 
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***** TSTR OUTPUT TEXT FILE ***** 
6 11 1987 
RUN NUMBER 1 


***** SOLUTION PARAMETERS ***** 

CC VOLUME = 1.30D-02 M**3 

CC HEAT TRANS. RATE => -6.50000D+07 J/S 

MIN ALLOWABLE MOLE NO. = 1.00000D-18 KGMOLE/KG 

INITIAL MASS OF LIQ IN THE COMB CHAMBER = 1.00000D-15 KG 

INITIAL DROPLET DIAMETER = 1.00000D-07 M 

MIN ALLOWABLE DROPLET DIAMETER = 5.000D-07 M 

EXIT TURBINE FLOW RATE CONSTANT = 9.233D-01 

HOT GAS MANIFOLD PRESSURE = 1.723D+06 PA 

NDEBUG (1 FOR PRINITING, 0 FOR SUPRESSING OUTPUT) 
1111111000 
0001000000 

DROPLET DIAMETER COEFFICIENT = 4.8500D-01 

AREA OF HYDROGEN INJECTOR ANULUS = 1.3290D-05 

AREA OF THE OXIDIZER INJECTOR = 1.0968D-05 

MCON (=0 FOR CONSTANT MASS FLOW; =1 FOR PIPING SYSTEM) 1 

KINET (=0 FOR NO KINETICS; =1 FOR CALL TO RADSS) 0 

INTEGRATION PARAMETERS: 

MITER, METH, NDGEAR 220 
IER, EPSI 10 5.00D-07 

INPUT STEP SIZE => l.OOD-12 SEC 
INITIAL TIME = O.OOD+OO SEC 
OUTPUT INTERVAL = 5.00D-0** SEC 

FINAL TIME = l.QOD-Ol SEC 
TIME SHUTDOWN BEGINS = O.OOD+OO 

FOR SS RATE APPROXIMATION: 

CONVERGENCE CRITERION (MAX % CHANGE) = 1.00D-06 

MAX NO ITERATIONS = 25 

ERROR CHECK IDGT = 0 

IMSL ERROR FLAG JER = 10 

OXIDIZER FEED SYSTEM INPUTS: 

NO. OF PIPES 11 NO. OF NODES 12 
NO. OF HELIUM ORIFICE PIPES 7 

PIPE IN WHICH THE INTERFACE IS INITIALLY LOCATED 8 
RFLIQ FOR THE OX PIPES = 1.303D-01 

RFGAS FOR THE OX PIPES = 1.827D-01 


RFLIQ FOR THE HE ORIFICE = 
RFGAS FOR THE HE ORIFICE = 
RFLIQ FOR THE HE VALVE 
RFGAS FOR THE Iffi VALVE 
INJECTOR AND MANIFOLD CD = 
UPSTREAM SOURCE PRESSURE = 
INITIAL BACKFLOW INTERFACE 


1.367D-01 
5.283D-01 
1.720D+04 
1.827D+03 
5.000D-01 
5.16964D+06 PA 
POSITION = O.OOOD+OO 


PIPE GEOMETRY: 


PIPE 

NODE 1 

NODE 2 

AREA (M**2) 

LENGTH (M) 

1 

1 

2 

8. 1073D-0S 

0.30000 

2 

2 

3 

8. 1073D-05 

0.30000 

3 

3 

4 

8. 1073D-05 

0.30000 

4 

4 

5 

8. 1073D-05 

0.30000 

5 

5 

6 

8.1073D-05 

0.30000 

6 

6 

7 

8.1073D-05 

0.30000 

7 

7 

8 

8. 1073D-05 

0.30000 

8 

8 

9 

2.0268D-03 

0.00500 

9 

9 

10 

2.0268D-03 

0.04000 

10 

10 

11 

2.0268D-03 

0.04000 

11 

11 

12 

2.0268D-03 

0.04000 


PIPE INITIAL CONDITIONS: 


PIPE 

VELOCITY (M/SEC) 

LIQUID/VAP I NT. POSITION (% LENGTH) 

1 

1.0000D-03 

l.OOOOD+OO 

2 

1.0000D-03 

1 . OOOOD+OO 

3 

1.0000D-03 

l.OOOOD+OO 

4 

1.0000D-03 

l.OOOOD+OO 

5 

1.0000D-03 

l.OOOOD+OO 

6 

1.0000D-03 

l.OOOOD+OO 

7 

1.0000D-03 

l.OOOOD+OO 

8 

1.0000D-03 

0 . OOOOD+OO 

9 

1.0000D-03 

0 . OOOOD+OO 

10 

1.0000D-03 

0. OOOOD+OO 

11 

1.0000D-03 

0. OOOOD+OO 

INITIAL NODAL DENSITIES: 


NODE 

DENSITY (KG/M**3) 


1 

2.0739D+01 


2 

2.0739D+01 


3 

2.0739D+01 


4 

2.0739D+01 


5 

2.0739D+01 


6 

2.0739D+01 


7 

2.0739D+01 


8 

2.07 39D+01 


9 

9.9070D+02 


10 

9.9070D+02 


11 

9.9070D+02 




» 

12 9.9070D+02 

***** REACTION rate data from npt ***** 


FOLLOWING REACTION IGNORED 

C3H8 02 > C2H4 H2 

REACTANT C3H8 NOT FOUND IN SPECIES LIST 

FOLLOWING REACTION IGNORED 

C8H18 02 > C2H4 H2 

REACTANT C8H1 NOT FOUND IN SPECIES LIST 

FOLLOWING REACTION IGNORED 

C2H4 02 > CO H2 

REACTANT C2H4 NOT FOUND IN SPECIES LIST 

FOLLOWING REACTION IGNORED 

C6H6 02 > CO H2 

REACTANT C6H6 NOT FOUND IN SPECIES LIST 

FOLLOWING REACTION IGNORED 

FUEL 0IL02 > CO H2 

REACTANT FUEL NOT FOUND IN SPECIES LIST 

FOLLOWING REACTION IGNORED 

H20(L) > H20 

REACTANT H20( NOT FOUND IN SPECIES LIST 

REACTION 1 02 H2 > OH OH 

DATA COMMENT: PRATT 

FWD RATE DATA: BX = 11.903 TEN = 0.000 TACT = 22661.000 

CALCULATED REVERSE RATE CONSTANTS 

STD DEVIATION = 2.428D+00 CORR COEFF = 9.475D-01 

REV RATE DATA : BX2 = 19.489 TEN2 = 0.000 TACT2 = 30793.276 

REACTION 2 OH H2 > H20 H 

DATA COMMENT: PRATT 

FWD RATE DATA: BX = 11.439 TEN = 0.000 TACT = 5187.000 

CALCULATED REVERSE RATE CONSTANTS 

STD DEVIATION = 4.078D-01 CORR COEFF = 9.939D-01 

REV RATE DATA : BX2 = 14.133 TEN2 = 0.000 TACT2 = 15784.008 

REACTION 3 02 H > OH 0 

DATA COMMENT: PRATT 

FWD RATE DATA: BX = 11.677 TEN = 0.000 TACT = 8712.000 

CALCULATED REVERSE RATE CONSTANTS 

STD DEVIATION = 1.892D+00 CORR COEFF = 8.664D-01 
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REV RATE DATA : BX2 = 17.639 TEN2 = 0.000 TACT2 = 14056.891 

REACTION 4 0 H2 > OH H 

DATA COMMENT: PRATT 

FWD RATE DATA: BX = 9.352 TEN = 0.000 TACT = 3903.000 

CALCULATED REVERSE RATE CONSTANTS 

STD DEVIATION = 5.361D-01 CORR COEFF = 9.459D-01 

REV RATE DATA : BX2 = 10.976 TEN2 = 0.000 TACT2 = 6690.385 

REACTION 5 0 H20 > OH OH 

DATA COMMENT: PRATT 

FWD RATE DATA: BX = 11.095 TEN = 0.000 TACT = 9115.000 

CALCULATED REVERSE RATE CONSTANTS 

STD DEVIATION = 1.283D-01 CORR COEFF = 9.217D-01 

REV RATE DATA : BX2 = 10.025 TEN2 = 0.000 TACT2 = 1305.377 

REACTION 6 H H M > H2 M 

DATA COMMENT: PRATT 

FWD RATE DATA: BX = 12.699 TEN = -1.150 TACT = 0.000 

CALCULATED REVERSE RATE CONSTANTS 

STD DEVIATION = 1.051D-01 CORR COEFF = l.OOOD+OO 

REV RATE DATA : BX2 = 11.955 TEN2 = 0.000 TACT2 = 51707.198 

REACTION 7 0 0 M > 02 M 

DATA COMMENT: PRATT 

FWD RATE DATA: BX = 9.672 TEN = -0.278 TACT = 0.000 

CALCULATED REVERSE RATE CONSTANTS 

STD DEVIATION = 1.274D+00 CORR COEFF = 9.941D-01 

REV RATE DATA : BX2 = 7.667 TEN2 = 0.000 IACT2 = 50190.523 

REACTION 8 0 H M > OH M 

DATA COMMENT: PRATT 

FWD RATE DATA: BX = 10.627 TEN = 0.000 TACT = -1400.000 

CALCULATED REVERSE RATE CONSTANTS 

STD DEVIATION = 6.110D-01 CORR COEFF = 9.988D-01 

REV RATE DATA : BX2 = 15.564 TEN2 = 0.000 TACT2 = 54467.238 

REACTION 9 H OH M > H20 M 


DATA COMMENT: PRATT 

FWD RATE DATA: BX = 10.778 TEN = 0.000 TACT = -252.000 

CALCULATED REVERSE RATE CONSTANTS 

STD DEVIATION = 4.828D-01 CORR COEFF = 9.995D-01 

REV RATE DATA : BX2 = 16.785 TEN2 = 0.000 IACT2 = 63424.861 

FOLLOWING REACTION IGNORED 

H 02 M > H02 M 

REACTANT H02 NOT FOUND IN SPECIES LIST 
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1 

2.1829D+01 

1 . OOOOD+OO 

2 

2.1829D+01 

1. OOOOD+OO 

3 

2. 1829D+01 

1. OOOOD+OO 

4 

2. 1829D+01 

1. OOOOD+OO 

5 

2. 1829D+01 

1. OOOOD+OO 

6 

2.1829D+01 

1. OOOOD+OO 

7 

2.1829D+01 

1. OOOOD+OO 

8 

8. 7429D-01 

1.0216D-01 

9 

8.7429D-01 

0. OOOOD+OO 

10 

8.7429D-01 

0. OOOOD+OO 

11 

8.7429D-01 

0. OOOOD+OO 


NODAL DENSITIES AND PRESSURES: 

NODE DENSITY (KG/M**3) PRESSURE (PA) 


1 

2.0739D+01 

5. 1696D+06 

2 

2.0253D+01 

5.0484D+06 

3 

1.9775D+01 

4.9292D+06 

4 

1.9306D+01 

4.8123D+06 

5 

1.8845D+01 

4.6976D+06 

6 

1.8395D+01 

4.5853D+06 

7 

1.795UD+01 

4.47S3D+06 

8 

1. 7522D+01 

4.3677D+06 

9 

9.9070D+02 

4.2120D+06 

10 

9.9070D+02 

4.2119D+06 

11 

9.9070D+02 

4.2119D+06 

12 

9.9070D+02 

4.2119D+06 
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5 

6 

7 

8 
9 

10 

11 

12 


1.8900D+01 
1.84720+01 
1. 8056D+O1 
1. 76510+01 
9.90700+02 
9.90700+02 
9.90700+02 
9.90700+02 


4.71120+06 

4.60450+06 

4.50070+06 

4.39990+06 

4.28550+06 

4.28550+06 

4.28550+06 

4.28550+06 


***** INTEGRATING ROUTINE CAT, LED ***** 


CALL FOR STEP 3 OUT OF 200 STEPS 

STEP BEGINS AT TIME = 1.00000D-03 AND ENDS AT TIME = 1.50000D-03 


***** CC AND INLET CONDTNS ON RETURN ***** 
***** FROM THE INTEGRATING SUBROUTINE ***** 


TIME FOLLOWING INTEGRATION = 1.50000D-03 SEC 

FINAL STEP SIZE = 1.75435D-05 SEC 

ERROR FLAG, IER = 0 


SPECIES 

NAME 

S2 (KGMOLE/KG) 

ISTRM 

SI (KGMOLE/KG) 

STATE 

1 

HE 

1.00944D-18 

3 

2.49838D-01 

G 

2 

H2 

4.82805D-01 

2 

4.96032D-01 

G 

3 

02 

1.03397D-11 

1 

3 . 12500D-02 

L 

4 

H20 

1.48005D-03 

0 

0.000000+00 

G 

5 

H 

1.00000D-18 

0 

0.000000+00 

G 

6 

0 

1.00000D-18 

0 

0.000000+00 

G 

7 

OH 

1.00000D-18 

0 

0.000000+00 

G 


DENSITY OF LIQUID IN THE CC = 
BACKFLOW INTERFACE POSITION = 
CC TEMPERATURE = 

CC PRESSURE 
CC DENSITY 
AVG. MOL. WT. = 

CC VOLUME 

OMV, FMV, PMV, EMV = 

UPSTREAM PRESSURE = 


5.07032D-03 KG 
0.000000+00 
465.0223 K 
4.210580+06 PA 
2.248730+00 KG/M**3 
4.84285D-01 KG/KGMOLE 
1.30000D-02 M**3 

1.3500+02 1. 615D+03 0.0000+00 

5.16960+06 PA 


INJECTOR PRESSURE DROP = 7.42060+02 PA 


1 . 694D+03 KG/M3/S 


STREAM 

TEMP (K) 

MDOT (KG/ SEC) 

SF (KGMOLE/SEC) 

SH (JOULES/SEC) 

1 

120.000 

1.755560+00 

0 . 0000000+00 

-1.1520830+06 

2 

160.000 

2.100000+01 

8.012821D+02 

-2. 5185900+09 

3 

120.000 

0.000000+00 

0.0000000+00 

0.0000000+00 

PIPE CONDITIONS: 





PIPE VELOCITY (M/SEC) LIQUID/VAP INT. POSITION (% LENGTH) 


SPECIES 

NAME 

S2 (KG’WLE/KG) 

1 

HE 

1.00027D-18 

2 

H2 

4.90960D-01 

3 

02 

1.51292D-11 

4 

H20 

5.67428D-04 

5 

H 

1.00000D-18 

6 

0 

1.00000D-18 

7 

OH 

1.00000D-18 


ISTRM SI (KGMOLE/KG) STATE 
3 2.49838D-01 G 

2 4.96032D-01 G 

X 3.12500D-02 L 

0 O.OOOOOD+OO G 

0 O.OOOOOD+OO G 

0 O.OOOOOD+OO G 

0 0 . OOOOOD+OO G 


DENSITY OF LIQUID IN THE CC = 
BACKFLOW INTERFACE POSITION = 
CC TEMPERATURE = 

CC PRESSURE 
CC DENSITY 
AVG. MOL. WT. = 

CC VOLUME 

OMV, FMV, PMV, EMV = 

UPSTREAM PRESSURE = 


1.11943D-03 KG 
O.OOOOOD+OO 
467.4076 K 
4.28538D+06 PA 
2 . 24344D+00 KG/M**3 
4.91528D-01 KG/KGMOLE 
1.30000D-02 M**3 

7.201D+01 1.615D+03 O.OOOD+OO 

5.1696D+06 PA 


INJECTOR PRESSURE DROP = 2.U02D+02 PA 


1.698D+03 KG/M3/S 


STREAM 

TEMP (K) 

MDOT (KG/SEC) 

SF ( KGMOLE/ SEC) 

SH (JOULES /SEC) 

1 

120.000 

9. 36176D-01 

0 . OOOOOOD+OO 

-6. 143656D+05 

2 

160.000 

2.10000D+01 

8.012821D+02 

-2.S18590D+09 

3 

120.000 

O.OOOOOD+OO 

0. OOOOOOD+OO 

0. OOOOOOD+OO 

PIPE CONDITIONS: 




PIPE 

VELOCITY (M/SEC) 

LIQUID/VAP INT. 

POSITION (% LENGTH) 

1 

1.1627D+01 

l.OOOOD+OO 



2 

1. 1627D+01 

l.OOOOD+OO 



3 

1. 1627D+01 

l.OOOOD+OO 



4 

1. 1627D+01 

l.OOOOD+OO 



5 

1.1627D+01 

l.OOOOD+OO 



6 

1. 1627D+01 

l.OOOOD+OO 



7 

1. 1627D+01 

l.OOOOD+OO 



8 

4.6623D-01 

3. 3539D-02 



9 

4.6623D-01 

O.OOOOD+OO 



10 

4.6623D-01 

O.OOOOD+OO 



11 

4. 6623D-01 

O.OOOOD+OO 



NODAL DENSITIES AND PRESSURES: 



NODE 

DENSITY (KG/M**3) 

PRESSURE (PA) 



1 

2.0739D+01 

5 . 1696D+06 



2 

2.0260D+01 

5.0501D+06 



3 

1.9793D+01 

4.9339D+06 



4 

1.9340D+01 

4.82 10D+06 




STREAM 

TEMP (K) 

MDOT (KG/SEC) 

SF (KGMOLE/SEC) 

SH (JOULES/SEC) 

1 

120.000 

2.71235D-01 

0.0000000+00 

-1.7799820+05 

2 

160.000 

2.100000+01 

8.0128210+02 

-2 . 518590D+09 

3 

120.000 

0.000000+00 

0.0000000+00 

0.0000000+00 

PIPE CONDITIONS: 




PIPE 

VELOCITY (M/SEC) 

LIQUID/ VAP INT. 

POSITION (% LENGTH) 

1 

3.3503D+00 

1.00000+00 



2 

3.35030+00 

1.00000+00 



3 

3 . 3503IH00 

1.00000+00 



4 

3 . 3503D+00 

1.00000+00 



5 

3.35030+00 

1.00000+00 



6 

3.35030+00 

1.00000+00 



7 

3.35030+00 

1.00000+00 



8 

1.3508D-01 

4.69740-03 



9 

1.3508D-01 

0.00000+00 



10 

1.3508D-01 

0.00000+00 



11 

1.3508D-01 

0.00000+00 



NODAL DENSITIES AND PRESSURES: 



NODE 

DENSITY (KG/M**3) PRESSURE (PA) 



1 

2.07390+01 

5.16960+06 



2 

2. 04270+01 

5.09180+06 



3 

2.01210+01 

5.01560+06 



4 

1.9822D+01 

4.94100+06 



5 

1. 95290+01 

4.86810+06 



6 

1.92430+01 

4.79670+06 



7 

1.89620+01 

4.72680+06 



8 

1.86880+01 

4.65830+06 



9 

9.90700+02 

4.59410+06 



10 

9.9070D+02 

4.59410+06 



11 

9.90700+02 

4.5941D+06 



12 

9.90700+02 

4.59410+06 




***** integrating routine called ***** 


CALL FOR STEP 2 OUT OF 200 STEPS 

STEP BEGINS AT TIME = 5.000000-04 AND ENDS AT TIME = 1.00000D-03 

***** CC AND INLET CONDTNS ON RETURN ***** 

***** FROM THE INTEGRATING SUBROUTINE ***** 

TIME FOLLOWING INTEGRATION = 1.00000D-03 SEC 

FINAL STEP SIZE = 1.024570-05 SEC 

ERROR FLAG, IER = 0 
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5 H 1.00000D-18 0 O.OOOOOD+OO G 

6 0 1.00000D-18 0 O.OOOOOIH-OO G 

7 OH 1.00000D-18 0 O.OOOOOD+OO G 

CC TEMPERATURE = 550.0000 K 

CC PRESSURE = 5.1696D+06 PA 

CC DENSITY = 2.27907D+00 KG/M**3 

AVG. MOL. WT. = 4.96031D-01 KG/KGMOLE 

CC VOLUME = 1.30000D-02 M**3 

STREAM NAME TEMP MDOT SF HS SH 

(K) (KG/SEC) (KGMOLE/SEC) (J/KGMOLE) (J/SEC) 

1 02 120.00 O.OOOOOD+OO O.OOOOOD+OO -2.73000D+05 O.OOOOOD+OO 

2 H2 160.00 2.10000D+01 8.01282D+02 -3.14320D+06 -2.51859D+09 

3 HE 120.00 O.OOOOOD+OO O.OOOOOD+OO -3.54160D+06 O.OOOOOD+OO 

\ 

***** INTEGRATING routine called ***** 

CALL FOR STEP 1 OUT OF 200 STEPS 

STEP BEGINS AT TIME = O.OOOOOD+OO AND ENDS AT TIME = 5.00000D-04 

***** CC AND INLET CONDTNS ON RETURN ***** 

***** FROM THE INTEGRATING SUBROUTINE ***** 

TIME FOLLOWING INTEGRATION = 5.00000D-04 SEC 

FINAL STEP SIZE = 3.17799D-06 SEC 

ERROR FLAG, IER — 0 

SPECIES NAME S2 (KGMOLE/KG) ISTRM SI (KGMOLE/KG) STATE 

1 HE 1.00146D-18 3 2.49838D-01 G 

2 H2 4.95236D-01 2 4.96032D-01 G 

3 02 2.18638D-11 1 3.12500D-02 L 

4 H20 8.90040D-05 0 O.OOOOOD+OO G 

5 H 1.00000D-18 0 O.OOOOOD+OO G 

6 0 1.00000D-18 0 O.OOOOOD+OO G 

7 OH 1.00000D-18 0 O.OOOOOD+OO G 

DENSITY OF LIQUID IN THE CC = 4.84393D-05 KG 

BACKFLOW INTERFACE POSITION = O.OOOOOD+OO 
CC TEMPERATURE = 493.6106 K 

CC PRESSURE = 4.59581D+06 PA 

CC DENSITY = 2.26077D+00 KG/M**3 

AVG. MOL. WT. = 4.95325D-01 KG/KGMOLE 

CC VOLUME = 1.30000D-02 M**3 

OMV, FMV, PMV, EMV = 2.086D+01 1.615D+03 O.OOOD+OO 1.681D+03 KG/M3/S 

UPSTREAM PRESSURE = 5.1696D+06 PA 

INJECTOR PRESSURE DROP = 1.7713D+01 PA 
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ft 


5 


6 


4 

7 


7 


6 


5 


5 

2 


0 


7 


6 


6 

3 


0 


8 


6 


5 

7 


0 


9 


5 


7 

4 


0 


ISIDE ARRAY 








RXN 

ISIDE(K,RXN) 






1 


0 -I 

-1 

0 0 

0 

2 



2 


0 -1 

0 

1 1 

0 ■ 

-1 



3 


0 0 

-1 

0 -1 

1 

1 



4 


0 -1 

0 

0 1 

-1 

1 



5 


0 0 

0 

-1 0 

-1 

2 



6 


0 1 

0 

0 -2 

0 

0 



7 


0 0 

1 

0 0 

-2 

0 



8 


0 0 

0 

0 -1 

-1 

1 



9 


0 0 

0 

1 -1 

0 ■ 

-1 



REACTION RATE DATA IN 

SI UNITS 




RXN 1 

MODE 



BX 


TEN 


TACT 

1 

1 

FWD 


11.903 


0.000 


22661.000 



REV 


19.489 


0.000 


30793.276 

2 

I 

FWD 


11.439 


0.000 


5187.000 



REV 


14.133 


0.000 


15784.008 

3 

1 

FWD 


11.677 


0.000 


8712.000 



REV 


17.639 


0.000 


14056.891 

4 

1 

FWD 


9.352 


0.000 


3903.000 



REV 


10.976 


0.000 


6690.385 

5 

1 

FWD 


11.095 


0.000 


9115.000 



REV 


10.025 


0.000 


1305.377 

6 

3 

FWD 


12.699 


-1.150 


0.000 



REV 


11.955 


0.000 


51707.198 

7 

3 

FWD 


9.672 


-0.278 


0.000 



REV 


7.667 


0.000 


50190.523 

a 

3 

FWD 


10.627 


0.000 


-1400.000 



REV 


15.564 


0.000 


54467.238 

9 

3 

FWD 


10.778 


0.000 


-252.000 



REV 


16.785 


0.000 


63424.861 

***** CC AND INLET INITIAL CONDITIONS 




SPECIE 

NAME S2 (KGMOLE/KG) 

ISTRM 

SI (KGMOLE/KG) STATE 

1 

HE 

1. 

00000D-18 

3 

2.49838D-01 G 

2 

H2 

4. 

96031D-01 

2 

4.96032D-01 G 

3 

02 

3. 

12500D-11 

1 

3. 12500D-02 L 

4 

H20 

5. 

55062D-11 

0 

O.OOOOOD+OO G 
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***** "Z" ARRAY OF THERMO PROPERTY POLYNOMIAL COEFFICIENTS ***** 


ft 


ft 



NAME 

COEFFICIENTS 





• 

HE 

2 . 50000D+00 
-7 .45375D+02 
O.OOOOOD+OO 

O.OOOOOD+OO 
9. 15349D-01 
O.OOOOOD+OO 

O.OOOOOD+OO 
2. 50000D+00 
-7.45375D+02 

O.OOOOOD+OO 

O.OOOOOD+OO 

9.15349D-01 

O.OOOOOD+OO 

O.OOOOOD+OO 

• 

H2 

3.10019D+00 
-6. 77380D+02 
5.52103D-09 

5. 11195D-04 
-1. 96294D+00 
-1.81227D-12 

5.26442D-08 
3.05 744D+00 
-9.88904D+02 

-3.49100D-11 

_2.67652D-03 

-2.29970D+00 

3.69453D-15 

-5.80991D-06 

• 

02 

3.62195D+00 

-1.20198D+03 

-6.76351D-09 

7.36183D-04 

3.61509D+00 

2.15560D-12 

-1.96522D-07 

3.62560D+00 

-1.04752D+03 

3.62016D-11 
-1.87822D-03 
4 . 30528D+00 

-2.89456D-15 

7.05545D-06 


H20 

2.71676 D+OO 
-2.99058D+04 
-2.96374D-09 

2.94514D-03 

6.63057D+00 

8.07021D-13 

-8.02244D-07 

4.07013D+00 

-3.02797D+O4 

1.02267D-10 

-1.10845D-03 

-3.22700D-01 

-4.84721D-15 
4. 15212D-06 

• 

H 

2.50000D+00 

2.54744D+04 

O.OOOOOD+OO 

O.OOOOOD+OO 

-4.59898D-01 

O.OOOOOD+OO 

O.OOOOOD+OO 
2 . 50000D+00 
2.54744D+04 

O.OOOOOD+OO 

O.OOOOOD+OO 

-4.59898D-01 

O.OOOOOD+OO 

O.OOOOOD+OO 

• 

0 

2.53430D+00 

2.92311D+04 

-3.26049D-09 

-1.24782D-05 

4.96286D+00 

1.01520D-12 

-1.25627D-08 

3.03094D+00 

2.91365D+04 

6.90299D-12 
-2.25259D-03 
2 . 60993D+00 

-6.37971D-16 

3.98245D-06 

• 

OH 

2.91312D+00 
3.96471D+03 
2.08436D- 10 

9 . 54182D-04 
5.42887D+00 
-2.33843D-13 

-1.90843D-07 
3.83655D+00 
3. 67158D+03 

1. 27308D-11 
-1.07020D-03 
4. 98055D-01 

2.48039D-16 

9.48498D-07 


***** REACTION RATE INDICIES AND DATA ***** 

® REACTION INDICIES 

RXN (J) ID(1,J) ID(2,J) ID(3,J) ID(4,J) 

1 3 2 7 7 

2 7 2 4 5 

3 3 5 7 6 

• 4 6 2 7 5 
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CH30 M > CH20 

REACTANT CH30 NOT FOUND IN SPECIES LIST 

FOLLOWING REACTION IGNORED 

CH20 M > CHO 

REACTANT CH20 NOT FOUND IN SPECIES LIST 


FOLLOWING REACTION IGNORED 

CH20 0 > CHO 

REACTANT CH20 NOT FOUND IN SPECIES LIST 


OH 


FOLLOWING REACTION IGNORED 

CH20 H > CHO 

REACTANT CH20 NOT FOUND IN SPECIES LIST 


H2 


FOLLOWING REACTION IGNORED 

CH20 OH > CHO 

REACTANT CH20 NOT FOUND IN SPECIES LIST 


H20 


FOLLOWING REACTION IGNORED 

CHO M > CO 

REACTANT CHO NOT FOUND IN SPECIES LIST 


FOLLOWING REACTION IGNORED 

CHO H > CO 

REACTANT CHO NOT FOUND IN SPECIES LIST 


H2 


FOLLOWING REACTION IGNORED 

CHO OH > CO 

REACTANT CHO NOT FOUND IN SPECIES LIST 


H20 


FOLLOWING REACTION IGNORED 

CHO 0 > CO 

REACTANT CHO NOT FOUND IN SPECIES LIST 


OH 


***** SPECIES PROPERTIES ***** 

MAJOR CONSTITUENTS: 

INDEX NAME MOLEC. WT. 

1 HE ^*.00260+00 

2 H2 2.0160IH00 

3 02 3.2000D+01 

U H20 1.8016D+01 


RATE DETERMINING RADICALS 
INDEX NAME MOLEC. WT. 

5 H 1.0080D+00 

6 0 1.6000D+0! 

7 OH 1.7008D+01 
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FOLLOWING REACTION IGNORED 

CO OH > 

REACTANT CO NOT FOUND IN SPECIES 

FOLLOWING REACTION IGNORED 

CO 0 M > 

REACTANT CO NOT FOUND IN SPECIES 

FOLLOWING REACTION IGNORED 

C02 0 > 

REACTANT C02 NOT FOUND IN SPECIES 

FOLLOWING REACTION IGNORED 

CH4 0 > 

REACTANT CH4 NOT FOUND IN SPECIES 

FOLLOWING REACTION IGNORED 

CH4 H > 

REACTANT CH4 NOT FOUND IN SPECIES 

FOLLOWING REACTION IGNORED 

CH4 OH > 

REACTANT CHA NOT POUND IN SPECIES 

FOLLOWING REACTION IGNORED 

CH3 CHO > 

REACTANT CH3 NOT POUND IN SPECIES 

FOLLOWING REACTION IGNORED 

CH3 CH20 > 

REACTANT CH3 NOT FOUND IN SPECIES 

FOLLOWING REACTION IGNORED 

CH3 0 > 

REACTANT CH3 NOT FOUND IN SPECIES 

FOLLOWING REACTION IGNORED 

CH3 OH > 

REACTANT CH3 NOT FOUND IN SPECIES 

FOLLOWING REACTION IGNORED 

CH3 H M > 

REACTANT CH3 NOT FOUND IN SPECIES 

FOLLOWING REACTION IGNORED 

CH3 02 > 

REACTANT CH3 NOT FOUND IN SPECIES 


H 

LIST 


C02 

LIST 


CO 

LIST 


CH3 

LIST 


CH3 

LIST 


CH3 

LIST 


CH4 

LIST 


CHO 

LIST 


CH20 

LIST 


CH20 

LIST 


CH4 

LIST 


CH30 

LIST 


FOLLOWING REACTION IGNORED 



OUTPUT FROM TSTR.PLT1 


TIME (MSEC) MOLE NUMBER (KGM0LES/KG*1.0I>+06) 


0.000 

l.OOD-15 

4.96D+02 

3.13D-08 

5.55D-08 

l.OOD-15 

l.OOD-15 

l.OOD-15 

0.50Q 

1.00D-15 

4.95D+02 

2.19D-08 

8.90D-02 

l.OOD-15 

l.OOD-15 

l.OOD-15 

1.000 

l.OOD-15 

4.91D+02 

1.51D-08 

5.67D-01 

l.OOD-15 

l.OOD-15 

l.OOD-15 

1.500 

1.01D-15 

4.83D+02 

1.03D-08 

1.48D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

2.000 

l.OOD-15 

4.72D+02 

7.01D-09 

2.65D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

2.500 

1.01D-15 

4.61D+02 

4.75D-09 

3.88D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

3.000 

1.01D-15 

4.52D+02 

3.24D-09 

4.98D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

3.500 

1.01D-15 

4.44EH-02 

2.23D-09 

5.86D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

4.000 

1.02D-15 

4.38D+02 

1.56D-09 

6.47D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

4.500 

1.01D-15 

4.36D+02 

1.11D-09 

6.77D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

5.000 

1.01D-15 

4.3SD+02 

7.98D-10 

6.78D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

5.500 

l.OOD-15 

4.38D+02 

5.83D-10 

6. 52D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

6.000 

1.01D-15 

4.42D+02 

4.31D-10 

6.03D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

6.500 

l.OOD-15 

4.48D+02 

3.21D- 10 

5.41D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

7.000 

l.OOD-15 

4.54D+02 

2.40D-10 

4.73D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

7.500 

l.OOD-15 

4. 59D+02 

1.79D-10 

4.11D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

8.000 

l.OOD-15 

4.64D+02 

1.33D-10 

3.63D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

8.500 

l.OOD-15 

4.66D+02 

9.83D-11 

3.34D+O0 

l.OOD-15 

l.OOD-15 

l.OOD-15 

9.000 

1.01D-15 

4.67D+02 

7.21D-11 

3. 270+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

9.500 

l.OOD-15 

4.66D+02 

5.25D-11 

3. 410+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

10.000 

l.OOD-15 

4.63D+02 

3.80D-11 

3.690+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

10.500 

1.01D-15 

4.60D+02 

2.74D-11 

4.050+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

11.000 

1.01D-15 

4.57D+02 

1.98D-11 

4.420+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

11.500 

1.01D-15 

4.54D+02 

1.43D-11 

4.740+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

12.000 

l.OOD-15 

4.52D+02 

1.03D-11 

4.980+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

12.500 

1.01D-15 

4.50D+02 

7.53D-12 

5.100+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

13.000 

1.01D-15 

4.50D+02 

5.50D-12 

5.110+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

13.500 

l.OOD-15 

4.51D+02 

4.03D-12 

5.010+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

14.000 

l.OOD-15 

4.53D+02 

2.97D-12 

4.840+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

14.500 

1.01D-15 

4.55D+02 

2.19D-12 

4.620+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

15.000 

l.OOD-15 

4.57D+02 

1.62D-12 

4.40D+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

15.500 

l.OOD-15 

'4.58D+02 

1.19D-12 

4. 200+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

16.000 

l.OOD-15 

4.60D+02 

8.76D-13 

4.070+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

16.500 

l.OOD-15 

4.60D+02 

6.43D-13 

4.000+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

17.000 

1.01D-15 

4.60D+02 

4.71D-13 

4.010+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

17.500 

l.OOD-15 

4.60D+02 

3.44D-13 

4.080+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

18.000 

1.01D-15 

4.59D+02 

2.50D-13 

4.190+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

18.500 

l.OOD-15 

4. 57D+02 

1. 820-13 

4.320+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

19.000 

l.OOD-15 

4.56D+02 

1.33D-13 

4.450+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

19.500 

l.OOD-15 

4.55D+02 

9.67D-14 

4.560+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

20.000 

l.OOD-15 

4.55D+02 

7.06D-14 

4.630+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

20.500 

1.01D-15 

4.54D+02 

5.16D-14 

4.660+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

21.000 

1.01D-15 

4.54D+02 

3.77D-14 

4.650+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

21.500 

l.OOD-15 

4.55D+02 

2.77D-14 

4.600+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

22.000 

1.01D-15 

4.56D+02 

2.03D-14 

4.530+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 

22.500 

1.01D-15 

4.56D+02 

1.49D-14 

4.450+00 

l.OOD-15 

l.OOD-15 

l.OOD-15 
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23.000 1.00D-15 

23.500 1.01D-15 

24.000 1.00D-15 

24.500 1.00D-15 

25.000 1.01D-15 


4.57D+02 1.09D-14 
4.57IH02 8.01D-15 
4.58D+02 5.87D-15 
4.58D+02 4.30D-15 
4.58D+02 3.14D-15 


4.38D+00 1.00D-15 
4.32D+00 1.00D-15 
4.28IHOO 1.00D-15 
4.27D+O0 1.00D-15 
4.28D+00 1.00D-15 


1.00D-15 1.00D-15 
1.00D-15 1.00D-15 
1.00D-15 1.00D-1S 
l.OOD-15 1.00D-15 
1.00D-X5 1.00D-15 
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OUTPUT FROM TSTR.PLT2 


ft 


ft 


TIME 

TEMPERATURE 

PRESSURE 

OX MDOT 

LIQ DEN 

DROP DIAM 

(MSEC) 

00 

(MPA) 

(KG/S) 

(KG/M3) 

(MICRONS) 

0.000 

550.00 

5.170 

0.0000 

7.69D-14 

0.500 

0.500 

493.61 

4.596 

0.2712 

4.84D-05 

1.914 

1.000 

467.41 

4.285 

0.9362 

1.12D-03 

6.382 

1.500 

465.02 

4.211 

1.7556 

5.07D-03 

11.863 

2.000 

477.70 

4.308 

2.4446 

1.15D-02 

16.710 

2.500 

497.47 

4.515 

2.9397 

1.79D-02 

20.582 

3.000 

518.24 

4.780 

3.1864 

2.19D-02 

22.964 

3.500 

535.86 

5.063 

3.2084 

2.26D-02 

23.795 

4.000 

547.72 

5.328 

3.0371 

2.02D-02 

23.114 

4.500 

552.48 

5.541 

2.7347 

1.59D-02 

21.207 

5.000 

549.80 

5.679 

2.3266 

1.09D-02 

18.262 

5.500 

540.22 

5.726 

1.8799 

6.44D-03 

14.813 

6.000 

525.28 

5.679 

1.4397 

3.38D-03 

11.296 

6.500 

507.36 

5.555 

1.0855 

1.68D-03 

8.423 

7.000 

489.16 

5.382 

0.8521 

9.15D-04 

6.508 

7.500 

473.24 

5.196 

0.7644 

6.99D-04 

5.736 

8.000 

461.60 

5.028 

0.8243 

8.46D-04 

6.086 

8.500 

455.41 

4.905 

1.0086 

1.40D-03 

7.354 

9.000 

454.85 

4.838 

1.2706 

2.47D-03 

9.202 

9.500 

459.14 

4.830 

1.5540 

4.06D-03 

11.246 

10.000 

466.86 

4.871 

1.8115 

5.93D-03 

13.166 

10.500 

476.31 

4.949 

2.0109 

7.62D-03 

14.734 

11.000 

485.81 

5.047 

2.1227 

8.72D-03 

15.707 

11.500 

493.91 

5.151 

2.1455 

8.98D-03 

16.039 

12.000 

499.56 

5.243 

2.0899 

8.44D-03 

15.758 

12.500 

502.19 

5.311 

1.9692 

7.35D-03 

14.947 

13-000 

501.76 

5.349 

1.8173 

6.04D-03 

13.840 

13.500 

498.68 

5.352 

1.6559 

4.80D-03 

12.615 

14.000 

493.71 

5.325 

1.5095 

3.82D-03 

11.470 

14.500 

487.85 

5.275 

1.3982 

3.17D-03 

10.574 

15.000 

482.09 

5.213 

1.3378 

2.83D-03 

10.058 

15.500 

477.31 

5.150 

1.3287 

2.79D-03 

9.929 

16.000 

474.12 

5.096 

1.3684 

2.99D-03 

10.172 

16.500 

472.83 

5.059 

1.4435 

3.42D-03 

10.692 

17.000 

473.38 

5.042 

1.5411 

4.00D-03 

11.395 

17.500 

475.46 

5.045 

1.6382 

4.66D-03 

12.117 

18.000 

478.54 

5.065 

1.7237 

5.27D-03 

12.774 

18.500 

482.02 

5.096 

1.7833 

5.73D-03 

13.257 

19.000 

485.32 

5.133 

1.8125 

5.97D-03 

13.522 

19.500 

487.95 

5.169 

1.8108 

5.95D-03 

13.556 

20.000 

489.61 

5.199 

1.7825 

5.73D-03 

13.383 

20.500 

490.15 

5.219 

1.7334 

5.37D-03 

13.041 

21.000 

489.64 

5.227 

1.6764 

4.94D-03 

12.621 

21.500 

488.27 

5.224 

1.6207 

4.54D-03 

12.197 


VAP RATE 
(KG/M3 S) 
0.00 
20.35 
66.55 
119.96 

165.78 
202.48 

225.56 
234.25 

228.56 

210.73 
182.44 

148.46 
113.41 

84.80 

65.86 

58.41 

62.20 

75.07 

93.38 

113.27 

131.76 
146.86 
156.37 

159.79 

157.29 
149.53 

138.74 

126.72 
115.43 

106.59 
101.52 

100.29 

102.73 

107.84 
114.72 

121.76 
128.15 

132.84 

135.46 

135.82 
134.17 

130.83 

126.75 

122.60 
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22.000 

486.36 

5.211 

1.5724 

4.22D-03 

11.819 

118.88 

22.500 

484.28 

5.191 

1.5397 

4.01D-03 

11.551 

116.24 

23.000 

482.36 

5.169 

1.5259 

3.92D-03 

11.423 

114.99 

23.500 

480.88 

5.148 

1.5302 

3.94D-03 

11.432 

115.07 

24.000 

480.02 

5.131 

1.5495 

4.07D-03 

11.558 

116.31 

24.500 

479.83 

5.120 

1.5800 

4.27D-03 

11.773 

118.41 

25.000 

480.25 

5.117 

1.6153 

4.50D-03 

12.032 

120.95 
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